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The  purpose  o£  this  research  was  to  continue  the 
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point  measurement  process  for  feedback.  The  effort  was 
motivated  in  part  by  the  concurrent  research  being  conducted 
at  the  Air  Force  Weapons  Laboratory/  Kirtland  AFB/  NM.  The 
research  was  conducted  with  the  primary  goal  of  developing  a 

neutral  particle  beam  controller/  but  the  issues  discussed 
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The  purpose  is  to  point  the  centroid  of  a  Neutral  Particle 
Beam  (NPB)  at  an  intended  target.  A  Multiple  Model  Adaptive 
Estimator  using  elemental  Meer  Filters  is  used  to  estimate  the 
centroid  of  a  NPB  model  as  a  one-dimensional  first-order 
Gauss-Markov  position  process.  The  MMAE  Meer  Filter  is  also 
used  to  estimate  the  beam  time  constant.  "Merge  Method"  of 
filter  pruning  is  used  to  limit  the  size  of  the  elemental  Meer 
filters.  A  bank  of  three  Kalman  filters  are  used  to  estimate 
the  states  of  the  target  which  has  a  variable  dynamics  driving 
noise  strength.  The  target  is  modelled  as  a  third-order 
Gauss-Markov  position  process.  A  Multiple  Model  Adaptive 
Controller  is  designed  using  LQG  methods,  and  true  states  are 
replaced  by  their  best  estimates  by  invoking  the  principle  of 
assumed  certainty  eguivalance.  MMAE  Meer  Filter  performance 
analysis  is  performed  for  an  uncontrolled  beam  and  for  a 
controlled  beam.  Controller  baselines  are  established. 


This  study  is  motivated  by  the  research  on  Neutral 
Particle  Beams  (NPB)  being  done  at  the  Air  Force  Weapons 
Laboratory  as  a  part  of  the  Strategic  Defense  Initiative. 
Most  previous  interest  in  the  Neutral  Particle  Beam  has  been 
its  use  as  a  weapon.  Figure  1.1  is  one  concept  of  a  neutral 
particle  beam  weapon.  Although  there  remains  a  strong 
interest  in  NPB  weapons,  the  Air  Force's  current  primary 
interest  in  the  NPB  is  to  use  it  to  determine  the  material 
components  of  a  target  (14].  This  is  to  be  accomplished 
through  the  nuclear  Interactions  of  the  particle  beam  with 
the  target.  As  it  is  necessary  to  direct  the  beam  onto  the 
target  to  perform  this  mission,  we  must  be  able  to  control 
the  beam.  Prior  to  being  able  to  control  the  beam,  we  must 
be  able  to  estimate  the  beam's  location.  Thus,  the  purpose 
of  this  study  is  to  continue  the  research  on  tracking  and 
control  of  the  NPB. 

One  approach  to  estimating  the  beam  location  involves 
injecting  laser  energy  into  the  beam.  Photons  absorbed  by 
the  beam  and  then  reemitted  are  collected  on  a  detector 
array.  The  pattern  of  the  photon  arrivals  is  used  to  infer 
information  about  the  beam's  location.  The  photon  events 
can  be  modelled  as  Poisson  space-time  point  processes. 

Snyder  and  Fishman  (291  developed  an  estimator  to 
incorporate  such  events  to  estimate  the  centroid  of  the 
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Figure  1.1  NPB  Concept  [20] 
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source  of  photon  events.  The  structure  of  the  Snyder 
Fishman  (SF)  filter  is  similar  to  that  of  the  Kalman  filter 
structure.  The  significant  difference  is  that  the  arrival 
times  of  events  in  the  SF  filter  are  not  known  a  priori. 
Santiago  [27]  showed  that  the  SF  filter  is  very  sensitive  to 
noise  induced  point  process  events  that  are  representative 
of  "dark  currents"  found  in  photon  detectors.  In  1962, 

David  Heer  in  his  doctoral  dissertation  [21]  developed  an 
estimator  that  filtered  out  the  noise  events.  The  Heer 
filter  is  based  on  a  bank  of  SF  filters.  Bach  SF  filter  in 
the  bank  assumes  a  different  sequence  of  signal-induced  or 
noise-induced  events.  Those  SF  filters  which  assume  the 
current  event  is  a  signal,  use  the  event  to  update  the 
filter's  estimate.  Those  SF  filters  which  assume  the 
current  event  is  a  noise,  ignore  the  event.  As  there  is  an 
ever-growing  number  of  sequences  of  signal  and  noise  events, 
this  requires  an  ever-growing  number  of  SF  filters.  Meer 
then  devised  a  method  to  limit  the  number  of  SF  filters 
through  the  "Best  Half"  method.  This  method  of  limiting  the 
number  of  SF  filters  was  replaced  by  the  "Merge"  method 
Incorporated  by  Zicker  in  1983  [32]. 

Once  the  beam  estimation  problem  had  been  solved,  the 
problem  of  beam  control  was  examined.  Zicker  conducted  a 
feasiblity  study  with  a  proportional  gain  controller  [32]. 
This  was  followed  by  the  design  of  proportional  plus 
Integral  (PI)  controllers  by  Moose  and  Jamerson  [8,231.  The 
type  one  control  provided  by  the  PI  controller  resulted  in 

3 


superior  tracking  performance.  Jamerson  also  Incorporated  a 
more  realistic  target  model  in  his  study.  This  target  model 


was  based  on  a  third  order  Gauss-Markov  position  process, 
rather  than  a  first  order  Gauss-Markov  position  process  as 
In  the  previous  studies.  PI  controller  performance  against 
the  new  target  was  found  to  be  sluggish.  Johnson  then 
designed  a  multiple  model  adaptive  controller  with  elemental 
PI  controllers  based  on  separate  assumptions  about  the 
strength  of  the  target  dynamics  driving  noise.  Each  of 
these  controllers  was  based  on  Linear  models  and  Quadratic 
cost  with  full-state  feedback  to  allow  LQ  synthesis 
techniques  (181  to  be  used  in  the  determination  of  constant 
controller  gains.  Then  by  applying  the  principle  of  assumed 
certainty  equivalence  [18],  the  full-state  feedback  was 
replaced  by  estimates  provided  by  the  Meer  beam  filter  and 
Kalman  target  filter. 

1.2  The  Problem 

In  examining  the  robustness  of  the  PI  controller  to 
variations  in  true  beam  parameters  (while  filter  assumed 
values  remained  constant),  Johnson  uncovered  a  significant 
robustness  problem  when  the  beam  time  constant  was 
mismodelled.  For  small  variations  in  the  true  beam  time 
constant,  the  Meer  filter  was  unable  to  estimate  the  beam 
position  precisely  and  the  controller  went  unstable.  Since 
stability  robustness  is  the  heart  of  any  feedback  controller 
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design,  this  problem  had  to  be  evaluated  and  a  solution 
found.  Johnson  evaluated  the  problem  and  one  of  his 
recommendations  was  to  apply  multiple  model  adaptive 
estimation  to  the  beam  time  constant  19). 


1 . 3  Scope 


This  document  assumes  the  reader  is  familiar  with  the 
basic  elements  of  control  theory  and  stochastic  processes. 
The  beam,  target  and  measurement  models  of  this  study  are 
limited  to  one  spatial  dimension.  This  is  done  to  allow 
simple  performance  characteristics  to  be  studied  without  the 
complications  of  cross  correlation  between  multiple  states 
and  parameters.  The  target  and  beam  are  assumed  to  be  well 
described  by  linear  models.  Where  reasonable,  random 
variables  or  distributions  are  assumed  to  be,  or  simplified 
to  be,  either  Gaussian  or  uniform  as  appropriate. 

Controller  design  is  based  on  LQ  synthesis  and  Multiple 
Model  Adaptive  Controller  techniques.  Implicit  Model 
Following  techniques  are  Investigated  as  a  method  to 
determine  controller  gains  which  offer  robustness  and 
performance.  Beam  actuators  are  assumed  to  be  linear 
devices  commanded  through  zero-order-hold  circuits  that 
maintain  the  commanded  feedback  until  the  following 
controller  update. 


il 
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This  project  begins  by  designing  a  multiple  model 
adaptive  estimator  using  elemental  Meer  filters  based  on 
separate  assumptions  of  the  beam  time  constant.  The 
weighted  state  estimates  are  combined  to  provide  an  adaptive 
state  estimate.  The  same  weightings  are  applied  to  the 
filter-assumed  beam  time  constants  to  provide  an  adaptive 
parameter  estimate.  The  weightings  are  based  on  the 
residual  performance  of  the  individual  filters.  The 
controller  uses  Multiple  Model  Adaptive  Controller 
techniques  to  design  an  adaptive  controller  based  on  a  bank 
of  PI  controllers.  Each  elemental  controller  is  based  on  a 
separate  combination  of  assumed  beam  time  constant  and 
strength  of  the  target  dynamics  driving  noise.  In  an 
alternative  approach,  the  parameter  estimate  from  the  MMAE 
Meer  is  used  to  propagate  the  state  estimate  of  a  separate 
Meer  filter  that  is  updated  with  the  adaptive  state 
estimate.  This  adaptive  estimate  is  then  supplied  to  a 
controller  based  on  a  nominal  beam  time  constant. 

As  in  previous  controller  designs,  the  constant 
controller  gains  are  determined  through  LQ  synthesis 
techniques  assuming  full-state  feedback.  However,  in  this 
study.  Implicit  Model  Following  is  used  as  an  aid  to 
determine  the  set  of  weighting  matrices  that  results  in  a 
controller  that  provides  robustness  and  performance.  The 
full-state  feedback  is  then  replaced  by  the  Meer  and  Kalman 


filter  estimates  through  the  application  of  assumed 
certainty  equivalence  (18]. 

The  estimator  and  controller  designs  are  validated  using 
computer  simulations.  Since  any  adaptive  filter  is 
nonlinear  and  the  resulting  filter-controller  combination  is 
also  nonlinear/  covariance  analysis  is  inappropriate  for 
this  simulation.  Instead/  Monte  Carlo  simulations  are  used 
to  collect  statistics  to  evaluate  performance.  Controller 
performance  is  measured  against  two  baselines.  The  first 
baseline  is  a  controller  with  full-state  feedback  and  access 
to  the  true  parameter  values.  The  second  baseline  is  a 
controller  receiving  estimates  from  a  single  beam  and  target 
filter  that  have  access  to  the  true  parameters.  The  first 
baseline  Indicates  the  best  possible  performance  and 
robustness  for  a  given  controller  design.  The  second 
baseline  provides  a  realistic  assessment  of  the  best 
performance  possible  from  controllers  that  operate  with 
state  estimates  rather  than  the  actual  states.  \ 

i 

< 

1.5  Summary  of  Remaining  Chapters 

Chapter  II  develops  the  estimator.  The  chapter  includes 
information  on  the  the  beam  model/  space  time  point 
processes.  Multiple  Model  Adaptive  Estimation/  the  Meer 
filter,  hypothesis  pruning,  and  the  MMAE  Meer  filter. 

Chapter  III  develops  the  controller.  Information  in  this 
chapter  includes  the  target  model,  LQ  synthesis.  Implicit 
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Model  Following,  and  various  controller  structures,  chapter 
IV  describes  the  performance  analysis  of  the  estimator  and 
controller.  This  includes  the  analyses  to  be  performed  and 
the  software  tools  used  to  generate  and  collect  the  data. 
Chapter  V  contains  the  results  of  the  analyses  on  the  MMAE 
Meer  filter  and  the  controllers.  Time  constraints  prevented 
the  analysis  of  a  controller  with  an  Implicit  Model 
Following  based  design.  Chapter  VI  contains  the  conclusions 
of  this  study  and  recommendations  for  subsequent  studies. 
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II .  Beam  Tracking 

The  first  step  In  tracking  and  pointing  of  a  Neutral 
Particle  Beam  (NPB)  Is  the  determination  of  the  beam's 
location,  or  tracking  the  beam.  In  his  doctoral  effort, 

Meer  [21]  developed  a  Multiple  Model  Adaptive  Estimator 
(MMAE)  using  elemental  Snyder-Fishman  filters  to  determine 
the  location  of  the  beam  centroid.  This  chapter  will  first 
describe  the  NPB  dynamics  model.  Next,  it  will  cover  the 
space-time  point  process  model  and  how  it  applies  to  the 
case  of  measurements  that  can  be  extracted  from  a  NPB  being 
Irradiated  with  laser  beams.  Then,  the  Snyder-Fishman 
filter  which  uses  that  space-time  point  process  as  its 
assumed  measurement  model  will  be  described.  The  Multiple 
Model  Adaptive  Estimator  will  be  presented  to  describe  the 
basic  structure  of  the  Meer  and  MMAE  Meer  filters.  Then, 
the  Meer  filter  will  be  described  which  uses  a  bank  of 
Snyder-Fishman  filters  to  estimate  the  NPB  centroid  in  the 
presence  of  noise  measurement  events.  Tree  pruning  and 
merging  will  then  be  described  as  methods  of  making  the  Meer 
filter  computationally  realizable.  Finally,  a  MMAE  using 
elemental  Meer  filters  will  be  developed  as  a  way  to 
accommodate  the  effects  of  an  uncertain  beam  time  constant. 

2 . 1  NPB  Model 

As  its  name  implies,  the  NPB  is  a  beam  of  particles.  In 
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modeling  the  NPB,  we  are  Interested  in  both  Its  nature  as 
Individual  particles  and  as  a  beam.  This  section  begins 
with  a  description  o£  the  particle  aspect  of  the  NPB.  This 
will  be  followed  by  a  description  of  the  beam  nature  that  is 
of  interest  to  us. 

The  term  neutral  particle  indicates  the  NPB  is  a 
collection  of  electrically  neutral  particles,  for  instance, 
monoatomlc  hydrogen,  which  is  composed  of  a  single 
negatively  charged  electron  orbiting  a  single  positively 
charged  proton.  This  gives  the  NPB  its  electrically  neutral 
characteristics  while  propagating  in  space  (201.  One 
approach  to  modelling  the  NPB  might  be  to  model  the  dynamics 
of  individual  particles  through  the  equation  [20:Chap  1,4-5) 

£  *  VmQ  (dv/dt  +  V2  v/c2  dv/dt  v)  (2-1) 

where 

f.  is  the  force  vector 

sIq  is  the  rest  mass  of  the  particle 

v  is  the  velocity  vector 

v  is  the  scalar  magnitude  of  the  velocity 

c  is  the  speed  of  light 

Y  is  a  relativistic  term  given  by 

Y  -  (  1  -  (v/c)2  1  ~1/2  (2-21 

Although  this  approach  might  be  highly  accurate,  to  model  a 

21 

lethal  NPB  would  require  on  the  order  of  10  such  models 
(25:59).  Fortunately,  McKee  indicates  that  for  the  portion 


of  the  MPB  that  we  are  Interested  In,  there  are  essentially 


no  collisions  between  particles  (20:8-11.  Since  there  are 


no  collisions,  a  group  of  these  particles  traveling  in 


essentially  the  same  direction  can  be  viewed  as  traveling 


along  a  long  tube  or  cylinder,  in  other  words,  viewed  as  a 


beam.  The  centerline  of  this  cylinder  describes  the 


centroid  of  the  beam.  The  centerline  also  describes  the 


direction  of  the  beam  velocity  vector. 


Since  the  objective  is  to  point  the  beam,  or  to  align 


the  velocity  vector  with  the  target,  the  movement  of  the 


beam  centroid  needs  to  be  modelled.  To  model  the  motion  of 


the  centroid,  we  must  consider  its  expected  response  to  our 


control  inputs  and  additionally  to  disturbance  inputs  and 


dynamics  driving  noises.  For  this  study  we  will  consider 


the  motion  of  the  beam  centroid  to  be  adequately  modeled  by 


the  linear  stochastic  differential  equation  (16:204): 


&£<t)  =  E.(t)2c(t)dt  +  &(t)&(t)dt  +  g(t)dfi(t)  ( 2-3a ) 


or,  in  the  less  rigorous  white  noise  notation  (16:204): 


d*(t)/dt  -  F(t)*(t)  +  fi(t)u(t)  +  £<t)w(t) 


( 2-3b) 


where 


&(t)  is  a  stochastic  vector  of  states  of  dimension  n 


characterizing  the  motion  of  the  beam  centroid. 


H.  &  is  the  beam  centroid  location,  as  indicated 


on  the  detector  array  (H.  &  is  a  projection  from 


Rn  onto  Rm;  m*2  for  a  plane,  m=l  for  this 
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experiment ) 

£L(t)  is  an  n-by-n  matrix  describing  the  system  dynamics 
fi.(t)  is  an  n-by-r  deterministic  input  matrix  relating 
the  control  inputs  to  the  beam  states 
y.(t)  is  an  r-vector  of  piecewise  constant  deterministic 
control  input  functions 

G(t)  is  an  n-by-s  noise  input  matrix  relating  the 
dynamic  driving  noises  to  the  beam  states 
&(t)  is  an  s-vector  Brownian  motion  with  statistics 
(E  is  the  expectation): 


Ejfi(t) 


E|tB(t)-B(t' ) Hfl(t)-a(t* ) r [  *  /  a(«)d« 


/:« 


( 2-4a ) 


with  Q.  an  s-by-s  matrix  of  piecewise  continuous 
functions  such  that  &(t)  is  symmetric  and  positive 
semldef inite . 

w(t)  is  a  s-vector  white  Gaussian  noise  (the  hypo¬ 
thetical  derivative  of  a  above)  with  statistics: 


Ejw(t) |  =  0 

E|w(t)w(t* )T\  -  a(t)«(t-t') 


( 2-4b) 


with  the  same  description  of  a  Just  given  (16:2041 


The  beam  centroid  states  are  propagated  using  Equation  (2-3) 


from  the  Initial  condition: 

2L<tQ)  *  JCg  (2-5) 

where  has  statistics  given  by: 

S|2L<t0)  }  ■  Sq  ( 2-6a) 

E|{x(t0)-m0](2t(t0)-a0)T|  -  E-o  (2-6b) 

where : 

tOq  is  the  mean  state  initial  condition 

Eg  is  an  n-by-n  covariance  matrix  that  is  symmetric  and 

positive  semidef inite 

For  pointing  of  the  NPB,  the  states  we  need  to  know  are 
the  location  and  the  direction  of  the  velocity  vector  or, 
alternatively,  the  centroid  location  at  two  points  along  the 
"cylinder"  at  essentially  the  same  time.  Either  case 
requires  one  to  gather  some  information  on  particle 
locations  in  the  beam.  One  approach  to  gathering  this 
information  is  to  shine  a  laser  beam  into  the  particle  beam 
to  "excite"  the  electrons  in  the  beam  into  higher  energy 
levels  or  quantum  states  (261.  When  the  beam  electrons 
"relax"  they  release  this  energy  in  the  form  of  photons  at 
specific  quantum  levels  (13:16-181  (See  Figure  2.1).  These 
photons  can  then  be  collected  by  using  a  photodetector 
array.  The  number  of  photons  that  arrive  in  a  given  amount 
of  time,  at  a  given  location  of  the  array  is  related  to  the 
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location  of  the  electrons  releasing  the  photons.  An 
electron  in  the  detector  that  is  energized  by  a  photon 
provides  a  signal  that  a  photon  event  has  occurred.  The 
model  for  the  location  and  time  of  arrival  of  these  photon 
events  is  covered  in  the  next  section. 


2 . 2  Beam  Measurement  Model 

This  section  models  the  signal  and  noise  processes 
occurring  at  the  detector  array.  Individual  signal  and 
noise  events  are  described  as  point  processes.  The  signal 
is  described  as  a  Poisson  space-time  point  process  with  an 
arrival  rate  density  that  has  a  Gaussian  spatial 
distribution.  Included  in  the  signal  model  is  an  array 
description  with  a  method  for  determining  appropriate  array 
size.  Then,  we  will  consider  the  presence  of  noise  at  the 
array.  The  noise  is  modeled  as  a  space-time  point  process 
uniformly  distributed  over  the  detector.  This  section 
concludes  with  expressions  for  combining  signal  and  noise 
process . 

2.2.1  The  Signal.  A  thorough  description  of  the 
emission  of  photons  from  a  NPB  and  their  detection  at  an 
array  is  not  necessary  for  this  study.  Therefore,  the 
objective  of  this  section  is  to  provide  only  enough 
information  to  justify  the  selected  model.  To  model  the 
signal,  we  would  like  to  know  the  rate  of  the  input  signal. 
Figure  2.2  depicts  the  emission  pattern  for  relativistic 
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Figure  2.2  Radiation  Intensity  versus  Emission  Direction  Tor  a  Unit  Charge  at  Various 
Velocities.  Forward  maximum  is  intensity  multiplication  relative  to  a  non-relativistic 
charge.  The  back  lobe  has  been  multiplied  by  the  given  factor  to  show  detail.  The  line 
labeled  with  an  angle  is  the  direction  of  xero  emission.  Q20) 


particles.  The  figure  Indicates  only  a  small  portion  of  the 
photons  emitted  In  the  NPB  will  arrive  at  a  detector  located 
to  the  side  of  the  beam.  As  described  in  Section  2.1,  the 
emission  of  photons  from  the  NPB  is  at  discrete  quantum 
levels.  These  same  quantum  effects  occur  at  the  detector 
surface.  As  a  photon  interacts  with  the  atoms  in  the 
detector,  an  electron  may  absorb  the  energy  of  the  photon, 
thereby  transitioning  to  a  higher  energy  level.  When  a 
transition  occurs,  an  event  is  declared  with  a  particular 
time  and  location.  Thus,  there  is  no  uncertainty  in  time  or 
location  of  the  event  occurrence.  Meer  and  Santiago 
[21,27,29]  further  described  the  arrival  of  individual 
photon  events  at  the  detector  array  as  a  Poisson  point 
process.  Meer  [21]  went  on  to  relate  this  signal  particle 
density  to  the  projection  of  the  beam  centroid  onto  the 
detector  array  as  (for  the  one  dimensional  case): 

i  (t,r,x(t))  «  r(t)exp( -Ir-H(t)x(t) ]TR”1(t)  [r-H(t )x( t ) 1/2 ) 

o 

(2-7) 

where 

X.  is  the  signal  rate:  the  arrival  of  photons  at  time  t 
at  position  r  on  the  detector,  given  that  the 
state  from  Equation  (2-3)  is  of  value  x(t)  so  that 
the  true  beam  centroid  is  at  H(t)x(t) 
r < t )  is  the  maximum  amplitude  of  the  rate  function 
r  is  the  spatial  location  of  the  signal  event  on  the 


detector 
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H ( t ) x ( t )  Is  the  signal-generated  beam  centroid  in  R1; 
a  projection  o£  the  beam  state  x(t)  onto  the  line 
of  detectors 

x(t)  is  a  state  stochastic  process  underlying  the  beam 
location  in  one  dimensional  space 
R  describes  the  spread  of  the  beam  and  is  therefore 
positive  definite.  The  square  root  of  R  is  the 
beam  halfwidth  (analogous  to  a  standard  deviation 
In  a  Gaussian  probability  density) 


Figure  2.3  depicts  the  arrival  rate  of  photons  or  signal 
particle  density  along  a  line  of  detectors  perpendicular  to 
the  beam  centerline.  As  can  be  seen  from  the  figure,  the 
tails  of  the  Gaussian  distribution  go  to  Infinity;  therefore 
only  an  infinitely  long  line  of  detectors  can  collect 
information  on  the  entire  density  function.  This  is 
physically  impossible  and  we  must  limit  the  length  of  the 
array.  If  a  is  the  standard  deviation,  cr  *%/R~ (the  beam 
half  width)  then  68.3%  of  the  photon  arrivals  occur  in  the 
Interval  (dropping  the  time  arguments)  [(Hx  -  c),(Hx  +  o)l. 
Similarly  99.7%  of  the  photons  arrivals  lie  in  the  interval 
( (Hx  -  3cr),(Hx  +  3cr)  ] .  Thus,  it  is  possible  to  define  a 
detector  size  to  collect  essentially  all  the  information  in 
the  density  function.  Santiago  (27:27-28)  defines  a  good 
"rule  of  thumb"  for  determining  detector  size  as  having  the 
detector  large  enough  such  that  the  distance  from  the 
centerpoint  of  the  Gaussian  rate  function  to  the  edge  of  the 
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detector  array,  is  at  least  six  beamwidths.  According  to 
Santiago  (27:28]  this  allows  the  beam  reasonably  large 
wanderings  while  still  allowing  edge  effects  to  be  ignored. 
This  results  in  12  beamwidths  or  24  beam  halfwidths  as  the 
appropriate  length  for  the  entire  line  of  detectors. 

Although  we  will  be  truncating  the  Gaussian  distribution,  we 
will  still  describe  the  the  arrival  rate  as  Gaussian  for 
mathematical  tractability.  In  reality,  the  array  will  also 
have  to  be  discretized  spatially  to  provide  position 
information.  This  will  limit  the  spatial  resolution  that  is 
possible.  However,  for  this  study,  the  array  will  be 
assumed  to  be  divided  into  infinitely  small  elements,  and 
thus  each  photoelectron  event  yields  both  a  time  and  exact 
location  of  the  event  on  the  detector  array.  Next  we  will 
consider  the  presence  of  noise  in  the  detector  and  model  it. 

2.2.2  The  Noise  Model.  The  sources  of  detector  noise 
are  basically  two-fold,  those  external  to  the  detector  and 
those  inherent  to  the  detector  itself  (13].  External  noise 
sources  will  be  covered  briefly,  followed  by  a  description 
of  inherent  noises.  A  total  noise  model  will  then  be 
described. 

External  or  background  noise  is  not  a  true  type  of 
noise,  but  refers  to  sensor  output  from  real  sources  of 
radiation  that  are  different  from  the  desired  source  (13). 
For  a  spaceborne  NPB,  some  of  the  background  noise  inputs 
are  the  sun,  stars,  earth  and  man-made  objects,  and  the 
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like.  Por  this  study  it  is  assumed  these  sources  are 
predictable  and  the  NPB  can  be  designed  to  eliminate  their 
impact,  or  that  their  effects  can  be  included  in  the  form  of 
models  to  be  used  for  inherent  noise  sources.  Therefore,  no 
separate  model  will  be  developed  for  the  external  noise 
sources . 

Unfortunately,  noise  sources  inherent  in  a  detector 
cannot  be  eliminated  and  therefore  must  be  modeled.  Types 
of  inherent  detector  noise  for  photodetectors  are 
Generation-Recombination  (G-R)  noise  and  Dark  Currents  [131. 
Here  we  have  assumed  that  only  detectors  capable  of 
detecting  individual  photons,  such  as  photomultiplier  tubes, 
avalanche  photodiodes,  or  semiconductor  photodiodes  will  be 
used  as  detectors.  G-R  noise  is  the  normal  generation  of  a 
free  electron  due  to  an  incident  photon  which  recombines 
within  the  semiconductor  lattice  prior  to  being  output  as  a 
signal  (131.  This  results  in  the  loss  of  a  signal  event. 

On  the  other  hand.  Dark  Current  is  the  output  of  an  apparent 
signal  event  that  was  not  generated  by  a  photon.  This  is 
due  to  the  statistical  distribution  of  electron  energies  in 
the  detector  materials  1131.  This  results  in  a  noise  event. 

Assuming  perfectly  uniform  detectors,  the  distribution 
of  the  noise  arrival  rate  will  also  be  uniform.  Therefore 
we  will  model  the  noise  arrival  rate  or  Poisson  point 
process  events  as  (21:1811: 


1A  =  constant  for  -L/2  s  r  s  L/2 
n 

(2-8) 

0  for  r  <  -L/2,  r  >  L/2 

where  An  is  the  noise  rate  and  L  is  the  length  of  the 
detector  array. 

Although  the  G-R  noise  is  clearly  dependent  on  the 
incoming  photons,  we  will  assume  the  effect  of  this  noise 
can  be  modeled  as  a  uniform  amplitude  reduction  in  the 
signal  arrival  rate  density  and  need  not  be  modeled 
separately.  As  the  noise  arrival  rate  due  to  Dark  Current 
is  independent  of  the  presence  of  photons,  we  will  assume 
the  noise  point  process  is  independent  of  the  signal  point 
process.  As  Meer  demonstrated,  the  signal  point  process  and 
noise  point  process  can  be  combined  into  a  single  Poisson 
point  process  with  rate  parameter  given  by  (21:21-231: 

A(t,r)  »  A_(t,r,x(t ) )  +  An(t,r)  (2-9) 

Often  in  signal  processing  that  involves  noise,  it  is 
desirable  to  consider  the  concept  of  a  signal  to  noise  ratio 
(SNR).  For  this  study  the  SNR  is  defined  as  the  average 
number  of  signal  induced  events  produced  for  every  noise 
event.  The  SNR,  thus,  can  be  expressed  as  (21): 

SNR  -  (p/(t£-t0) l/(An(t,r)  L)  ( 2-10a ) 


/2kR  r  /(An(t,r)  LI 


(2-10b) 


where 


p  is  the  number  of  signal-induced  events  occurring 
over  the  period  t£  -  tg 

An(t,r)  is  the  noise  arrival  rate  per  unit  length  of 
the  detector  array  given  in  Equation  (2-10) 

L  is  the  length  of  the  detector  array 

Figure  2.4  gives  a  graphical  comparison  of  the  signal  and 
noise  models  at  a  selected  signal  to  noise  ratio  upon  an 
appropriately  sized  detector.  The  next  step  is  to  develop 
an  efficient  method  to  recover  the  beam  location  Information 
from  the  signal  and  noise  information  using  the  models  we 
have  developed. 

2.3  Snvder-Flshman  Filter 

To  determine  the  centerline  of  the  beam  on  the  array  we 
need  to  use  an  estimator  that  can  optimally  process  all  the 
photon  measurements  that  we  provide  it.  It  must  be  able  to 
take  full  advantage  of  the  Gaussian  model  we  are  assuming 
for  spatial  distribution  of  signal-induced  photon  arrival 
rate  at  the  detector.  In  addition,  it  must  be  able  to 
process  events  that  arrive  at  random  intervals.  In  1975 
Snyder  and  Fishman  described  a  filter  to  "estimate  the 
centroid  of  a  swarm  of  fireflies  based  on  their  light 
flashes."  (29).  In  their  work,  they  described  the  firefly 
flashes  as  Poisson  point  processes.  Based  on  this,  they 
developed  a  minimum  mean  square  error  estimator  to  determine 
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Figure  2.4  Signal  and  Noise  Rate  Functions 


the  centroid  of  the  swarm.  The  estimator  Is  described  by 
the  following  differential  equations  (9:14-15,29): 


d&(t)  *  £(t)2t(t)dt  +  fi(t)u(t)dt 

)l£.-ii(t)x(t))  N(dt  x  dt)  (2-11) 

'R^ 


+  j  K(t: 


dE.(t )  =  F(t)P(t)dt  ♦  P(t)FT(t)dt  ♦  £(t)£T(t)dt 


I 


K(t)H.(t)E(t)N(dt  x  dj.)  (2-12) 


m 


K  ( t )  =  E.(t)HT(t)(H(t)E(t)|iT(t)  *■  b.(  t )  1_1  (2-13) 

where 


x(tQ)  =2^  and  P(tQ)  =  P^  are  the  Initial  conditions 
P(t)  is  the  filter-computed  error  covariance 
£(t)GT(t)  *  G(t)ft(t)£T(t);  with  £(t)*t  (See 
Equation  (2-4)  ) 

K(t)  is  the  filter  gain 

The  notation  found  In  Equations  (2-11)  and  (2-12)  Involves 
counting  Integrals,  where: 


f(t,t)H(dt  x  dt) 


/ 


0 


E  f(t1,t1)  «(t,t4) 

1*1 


Nfc  <  0 

(2-14) 


Nt  i.  1 


and  6(t,t4)  Is  a  Kronecker  delta.  Simply  stated,  if  no 
events  are  detected,  then  the  integral  equals  zero.  If  an 
event  Is  detected,  the  Integral  causes  a  jump  discontinuity. 
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t(t,£.)#  to  be  added  to  the  solution  of  the  differential 
equation  for  time  tj  (291. 

The  filter's  implementable  propagation  and  update 
equations  can  be  derived  from  Equations  (2-11)  and  (2-12). 
The  propagating  equations  for  use  between  (signal-induced) 
events  are: 

dx(t)  =  £(t)£.(t)dt  +  B(t)u(t)dt 
dP(t )  =  £(t)E(t)dt  +  P(t)£T(t)dt  +  G(t)GT(t)dt 
or  in  the  form  of  stochastic  difference  equations: 


x(ti+l}  *  + 


(2-17) 


E(tj  +  i)  =  ft(t1+1,t1)E.(tJ)  ft  (tl+1,t4)  +  SW  (2-18) 
where : 

— (ti  +  l'ti)  *s  the  state  transition  matrix  associated 
with  £ 

Qg  is  the  equivalent  discrete-time  noise  strength 
calculated  by: 


Bd(t 


1+1 


£(ti+l/«)£(a)a(a)GT(a)iT(ti+1/«)da  (2-19) 


where  ft  ■  £.  When  an  event  has  been  detected,  a  measurement 
update  can  take  place  as  defined  by  the  following  equations: 


$ 

(2-15) 

j£u 

(2-16) 

ft 

2L(tJ)  -  &<t")  +  jCttj  JlLj-ftCtj  )2L(t")  ] 


E(t*)  -  P(tj)  -  |C(ti)H(t1)E(t1)l 


where  the  Snyder -Fishman  gain  &  is  defined  by: 


(2-20) 


(2-21) 


K ( 1 1 )  «  E(t^)ttT(t1)(ii(t1)E(t^)HT(t1)  ♦  E(tj )  l'1  (2-22) 

For  this  study  we  will  use  a  one  dimensional  beam  model. 
Thus,  all  the  previous  equations  are  reduced  to  the  scalar 
case.  For  the  beam  dynamics  model  (F)  in  Equation  (2-3),  we 
will  assume  a  first-order  Gauss-Markov  position  process,  and 
thus : 


P  -  -1/7%, 


A  first-order  Gauss-Markov  process  often  provides  an 


(2-23) 


adequate  description  of  bandlimited  processes.  Figure  2.5 
depicts  the  nature  of  a  first-order  Gauss-Markov  process. 

The  equations  for  this  filter  look  remarkably  similar  to 
the  Kalman  Filter  equations  (19,21,27,29).  The  primary 
difference  is  that,  the  time  that  measurements  will  arrive 
is  not  known  a  priori.  In  1978,  Santiago  used  the 
Snyder-Fishman  filter  to  investigate  the  performance  of 
optical  trackers  when  operating  at  low  signal  rates  (27). 
Santiago  found  that  the  Snyder-Fishman  filter  performance 
deteriorated  markedly  in  the  presence  of  noise.  The  noise 
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Figure  2.5  First  Order  Gauss-Markov  Process 


in  this  case  was  that  due  to  dark  currents  within  the 


photodetector  as  well  as  background  noise.  In  an  ef£ort  to 
enhance  the  Snyder-Fishman  filter  performance,  Santiago  1271 
investigated  the  use  of  residual  monitoring  techniques  to 
perform  preliminary  rejection  of  noise  events.  In  an 
alternate  approach,  Meer  [21]  developed  the  Meer  filter  by 
applying  multiple  model  adaptive  estimation  techniques  to 
this  problem.  The  next  section  will  describe  multiple  model 
adaptive  estimation. 

2 . 4  Multiple  Model  Adaptive  Estimation 

As  mentioned  in  the  last  section,  Meer  applied  Multiple 
Model  Adaptive  Estimation  (MMAE)  techniques  to  solving  the 
problem  of  estimating  the  beam  position  with  Snyder-Fishman 
filters  in  the  presence  of  noise,  adaptively  determining 
whether  each  photoelectron  event  is  signal-induced  or  noise- 
induced.  MMAE  techniques  are  applied  to  the  problem  of 
estimating  the  beam  centroid  location  with  an  uncertain  beam 
time  constant,  adaptively  determining  the  true  value  of  the 
beam  time  constant.  This  section  will  present  a  development 
of  the  MMAE  prior  to  presenting  the  applications  of  this 
technique  in  later  sections. 

One  method  of  developing  a  MMAE  of  the  Meer  filter  is  to 
use  a  Bayesian  approach.  Following  the  development  by 
Maybeck  [17:129-136],  the  Bayesian  state  estimate  can  be 
expressed  as  the  expected  value  of  the  state,  conditioned  on 


the  measurement  history  2.^: 

X(tJ)  =  Elxtt^J^l  =  /-£*x(t)|1(£|l1)d£  (2-24) 

where  £  is  a  dummy  variable  of  integration  of  &(t),  and 
£x(t)  JZ(£l^i)  £s  the  con<**tlonal  probability  density 
function  of  &(t),  given  that  the  realization  2.^  of  the 
measurement  history  2  has  been  observed  up  to  the  current 
time.  If  we  let  a(t)  denote  the  vector  of  uncertain 
parameters  (beam  time  constant  for  example),  and  assume  §.(t) 
can  take  on  any  value  in  the  continuous  range  of  A,  the  last 
equation  can  be  rewritten  to  include  a(t): 

<2-«> 

where  a  is  the  dummy  variable  used  to  integrate  a(t)  over 
the  range  of  A.  According  to  Bayes'  rule,  the  conditional 
probability  density  function  can  be  expressed  as: 

£x<t>,„|Z(a'aIV  *  £x(t>  |a,Z  <ala'V  £a  |  Z  <aIV  ,2-2S' 

The  first  density  function  on  the  right  hand  side  of  the 
last  equation  is  the  density  of  the  state  (equal  to  the  beam 
centroid  location  since  H  *  1) ,  given  the  beam  time 
constant  and  measurement  histories.  It  has  a  mean  of 
j£(  t ^  )  and  a  covariance  of  E(t^)  as  computed  by  a  Meer 
filter  based  upon  a  particular  £  or  beam  time  constant.  The 
second  term  is  the  conditional  probability  of  £  conditioned 


on  the  measurement  history.  According  to  Bayes'  rule,  it 
can  be  expressed  as: 


£a|Z(-l2-i>  =  £a|z(i),Z(i-l)  ^l^i'^i-l* 

£a,z(i) |Z(i-l) <fiL l^i'^-i-l * 

s 

£Z(i)|Z(i-l)(|iil-i-l) 

£z(i)  |a,Z(i-l)  ^i  l-,-i-l)  £a|Z<-l2i-l> 

JA£z(i)|a,Z(i-l)(liil-/2-i-l)  £a|Z(-l-i-l,d“2_27) 


where  £ 


of 


Ez(  i)  |a,Z(  i-1)  ^i  I— '^1-1^  is  Gaussian'  with  a 
ft(t)2c(t“)  and  a  covariance  of  I H ( t ) EL< t ^ ) ftT < t ^ )  <■  RCtj)], 

again  as  produced  by  a  Heer  filter  based  on  a  specific  value 

of  &.  Here  £(i)  is  the  current  measurement  and  Z^^  is  the 

history  of  all  previous  measurements. 

The  last  equation  could  conceptually  be  evaluated 

iteratively  starting  from  £  (a).  By  combining  equation 

(2-25)  with  equation  (2-26)  and  switching  the  order  of 

integration,  we  can  generate  the  state  estimate  from: 


-l.'L 


M  £x(t)|a,Z(*l*'Vd£  1  £a|Z«MVd*  <2'28) 


The  term  within  the  brackets  of  Equation  (2-26)  is  the 
output  from  a  single  filter.  Unfortunately,  the  outer 
Integration  implies  an  uncountably  infinite  number  of  such 


filters  (one  for  every  possible  beam  time  constant).  This 
will  make  this  approach  computationally  impractical,  so  we 
let  the  uncertain  parameter  vector  (again  in  this  case  the 
beam  time  constant)  assume  only  a  finite  number  of  values. 


such 


as  the  discrete  vector  set  [g.^,  &2,  ...  ^J,  where 


is  the  number  of  elemental  filters  in  the  MMAE.  Thus,  the 
integration  is  replaced  by  a  finite  summation.  This  becomes 
the  foundation  of  the  multiple  model  adaptive  structure, 
where  each  estimate  within  the  bank  of  estimators  is  based 
on  an  assumed  discrete  parameter  value,  in  this  case  a  beam 
time  constant. 

To  adopt  the  Bayesian  approach  to  a  multiple  model 
structure,  an  a  priori  density  function  is  required  for 
^(tg)  (to  serve  as  an  initial  condition  on  the  Iteration 
to  be  discussed  next): 


£a(t0)(^  *  E  Pk(V 
U  k=l 


(2-29) 


where  p^ttg)  Is  the  probability  that  &  assumes  the  value 
at  time  tQ.  Therefore,  the  hypothesis  conditional  prob¬ 
ability,  p^t^),  is  a  recursive  relationship  expressed  as: 

pk(ti>  *  proMa.*^  |Z.(tl 


£z(l)  ja,Z(i-l) )  (141  lftk,Z-i— 1}  Pk^i-l* 
:j«l  £z(i)|a,Z(i-l)(liilaj'Z-i-l)  Pj^i-l1 


(2-30) 
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where  £z(i)  ja#Z(  i-1)  can  be  evaluated  as: 


Ez(i)  Ja,Z(i-l) 


I(2n)m*/2|Akj1/2j'1  expl-(l/2)£j(ti)&k(ti)~1tk(ti)  1 

(2-31) 


where 

r^(t1)  Is  the  residual  in  the  k-th  filter; 


tk(ti)  *  5-i  '  ttk(ti)ik(tI)  (2-32) 

^k(ti)  =  Hk(ti,Ek(ti)ak(ti)  +  Bk(tl)  (2-33) 

m  is  the  number  of  measurements  at  time  t,;  the 
dimension  of  ^(t^. 

Therefore,  the  conditional  mean  is  the  sum  of  the 
probabilistically  weighted  state  estimates  from  each 
estimator  within  the  bank: 


£(t*)  -  BCjc(tl )  |Z.(tl  )-Z.i  1 


( 2-34a) 
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Figure  2.6  depicts  the  structure 
The  preceding  discussion  was 


of  the  MMAE  estimator, 
presented  in  terms  of  a 
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multiple  model  structure  of  elemental  Meer  filters  to 
address  adaptation  on  beam  time  contant  rfi.  Applying  this 
same  structure,  Heer  developed  an  adaptive  filter,  using 
elemental  Snyder -Fishman  filters,  to  estimate  the  NPB 
centroid  despite  the  presence  of  noise  in  the  detector. 

2.5  Heer  Filter 

In  1982,  David  Meer  developed  an  adaptive  estimator  to 
solve  the  problem  of  estimating  the  NPB  in  the  presence  of 
noise  (21).  This  section  will  begin  by  presenting  the  basic 
Meer  filter.  Then,  the  methods  applied  by  Meer,  and  later 
Zicker,  to  reduce  the  computational  load  of  a  simple  Meer 
filter,  will  be  presented. 

2.5.1  Basic  Meer  Filter.  The  basic  Meer  filter  Is  an 
application  of  the  MMAE  structure  with  elemental 
Snyder -Fishman  filters.  Each  model  (Snyder-Fishman  filter) 
in  this  structure  is  based  on  a  separate  hypothesis  sequence 
for  the  order  of  incoming  signal  and  noise  events.  The 
order  of  the  occurrence  of  signal  and  noise  events  specifies 
each  individual  model's  hypothesis  sequence  as  shown  in 
Table  I.  For  each  event  labeled  as  a  signal,  the 
information  (location  of  the  event  on  the  detector  array  at 
that  time)  is  used  to  update  the  beam  location  estimate  of 
the  filter  based  on  that  hypothesis  sequence.  This  update 
is  performed  using  Equations  (2-20)  and  (2-21).  For  each 
event  labeled  as  a  noise,  no  update  occurs  so  that: 


Table  2.1  Hypothesis  Sequences 


Number  of  Events 


Possible 

Sequences 


A 

where  is  the  conditionally  expected  value  of  the  beam 
states  for  the  j-th  hypothesis  and  ^  is  the  variance  for 
the  j-th  hypothesis  sequence.  The  state  estimate  of  each 
elemental  filter  is  expressed  as: 

Xj(t)  =  Elx(t)  |h!jt,zNt]  (2-37) 

where  2Lj(t)  is  the  expected  value  of  the  beam  states,  x(t), 
conditioned  on  the  j-th  assumed  hypothesis  time  history 
h^fc  and  the  observations  history,  ZNt,  of  events  (t^,i^), 
where  i  *  l,2,...,Nfc. 

The  sequences  or  hypothesis  time  histories  can  be  formed 
into  conditional  probability  trees  to  illustrate  Bayes*  rule 
(See  Figure  2.7).  This  study  will  use  the  convention  that 
even  numbered  hypothesis  sequences  have  assumed  the  most 
recent  event,  or  M^-th  event,  was  a  signal.  For  a  hypo¬ 
thesis  sequence  corresponding  to  the  assumption  that  the 
most  recent  event  is  signal-induced,  the  conditional 
probability  is  given  by  (j  ■  2,4,...2Nt): 


Pr(h4|Z~)  =  Pr (s )Pr (h2 |Z  ) 


«  Pr (n)Pr (hjJZ1) 


=  Pr(s)Pr(h1|Z1) 


=  Pr(n)Pr(h1|Z1) 


Figure  2.7  Conditional  Probability  Tree 


where  there  are  now  2j  signal  related  hypothesis  sequences 
as  shown  in  Figure  2.7  and 
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(2-39) 


^8(t»t'INt'X(tNt,)  *  ^•n(tNt/tNt) 


where  A  and  A  are  the  signal  and  noise  rate  parameters 
s  n 


previously  described  in  Section  2.2.  The  weighting 
probability  starts  at  tQ/  with  the  initial  conditional 
probability  defined  as  Pr(hQ|Z°)  =  1. 

For  hypothesis  sequences  assuming  the  Nfth  event  was 
noise  induced,  the  conditional  probability  is  given  by 
(j  =  l,3,...2Nt-l): 


Pr(h2j-llzNt)  *  Pr(n)  p'<hj|zNt_1> 


(2-40) 


where 


Pr  (n) 


^n(tNfc'rNt  ] 


(2-41) 


^s(tHt/rNt'X(tNt)  }  +  <>ln(tNt'rNt) 


The  overall  state  estimate  of  the  MMAE  is  the  probabi¬ 
listically  weighted  average  of  individual  filter  state 
estimates  given  by: 


2Nt 


i(t)  -  Pr(h^fcjZNt)  < t ) 


(2-42) 
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with  an  MMAE  computed  error  variance  de£ined  as: 


L<t) 


-  E  Prth^lZ1**)  <P,(t)  +  (x  .(t)-jt(t)  1  lx  .(t)-x(t)  ]T) 
j«l  3  '  33  3  (2-43) 


The  structure  of  the  Meer  Filter  is  depicted  in  Figure  2.8. 

Note  that  the  number  of  hypothesis  sequences,  and 
therefore  the  number  of  models  or  Snyder-Fishman  filters, 
required  to  estimate  the  beam  location  optimally  by  the 
occurrence  of  events  is  2Wt.  This  quickly  becomes  a 
computationally  burdensome,  and  in  fact  intractable,  filter. 
For  example,  ten  events  requires  21®  *  1024  models  to 
implement  the  basic  filter  fully.  Meer  recognized  this 
drawback  of  the  basic  filter  design  (21)  and  proceeded  to 
develop  a  method  to  limit  the  expansion  of  the  hypothesis 
tree  by  pruning  some  of  the  hypothesis  sequences. 

2.5.2  Pruning  the  Hypothesis  Tree.  A  way  to  limit  the 
growth  of,  or  prune,  the  hypothesis  tree  was  needed.  Meer 
proposed  a  number  of  possible  methods  to  prune  the 
hypothesis  tree  and  chose  the  "Best  Half  Method"  as  the  most 
viable  (21).  Zicker  (32),  applied  the  work  of  Weiss, 
Upadhyay,  and  Tenney  (30,31)  to  limit  the  expansion  of  the 
hypothesis  tree  by  combining  the  branches  using  the  "Merge 
Method" . 

2. 5. 2.1  Best  Half  Method.  Meer  limited  the 
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expansion  of  the  hypothesis  tree  by  limiting  the  number  of 
events,  or  the  history  of  events,  that  each  hypothesis 
sequence  is  allowed  to  consider.  This  number  of  events 
defines  the  memory  depth  D  of  the  suboptimal  Heer  filter. 

For  this  study,  the  following  convention  is  used  to  describe 
the  age  of  an  event. 

When  an  event  A  first  occurs  it  has  an  age  of  1. 

When  the  next  event  occurs  Event  A's  age  becomes 
2.  Event  A  will  continue  to  age,  as  new  events 
come  in,  until  its  age  equals  the  depth.  When 
the  age  of  Event  A  reachs  depth,  it  is  the 
oldest  event.  When  the  next  event  occurs.  Event 
A  is  discarded. 

once  the  Meer  filter  reachs  the  prescribed  depth,  a 
comparison  is  made  between  the  probability  of  the  oldest 
event  being  a  signal  and  the  probability  the  oldest  event 
being  a  noise.  Conceptually,  if  the  comparison  indicates 
the  oldest  event  was  most  probably  a  signal,  then  the 
sequences  associated  with  the  oldest  event  being  a  signal 
(the  upper  half  of  the  tree)  are  retained  (See  Figure  2.9a). 
If  the  comparison  indicates  the  oldest  event  was  most 
probably  a  noise,  then  the  sequences  associated  with  the 
oldest  event  being  a  noise  (the  lower  half  of  the  tree)  are 
retained  (See  Figure  2.9b).  The  oldest  event  and  its 
probabilities  are  then  replaced  by  the  event  and 
probabilities  of  the  next-oldest  event.  The  process  of 
selecting  the  best  half  would  then  be  repeated  as  each  event 
comes  in.  Thus,  only  hypothesis  sequences  associated  with 
the  most  recent  "depth  D"  number  of  events  are  retained. 
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Hypothesis  sequence 
prior  to  pruning 


Hypothesis  Sequence 
after  pruning 


Figure  2.9a  Best  Half  Pruning  with  Signal  Most 
Probable.  (Depth  *  3) 


Hypothesis  sequence 
prior  to  pruning 


Hypothesis  Sequence 
after  pruning 


Figure  2.9b  Best  Half  Pruning  with  Noise  Most 
Probable.  (Depth  =  3) 


Although  the  complete  hypothesis  sequences  are  lost,  the 
portion  o£  the  sequence  that  Is  discarded  still  has  an 
Influence  on  the  hypothesis  time  history  (l.e.  which  set  of 
hypothesis  sequences  are  retained).  Conceptually,  although 
an  event  only  directly  impacts  the  estimation  of  the  beam 
state  for  "depth  D"  events  (while  it  is  in  the  filter),  it 
indirectly  impacts  state  estimation  after  leaving  the  filter 
by  determining  which  sequence  will  be  retained  for  future 
beam  state  estimation. 

2. 5. 2. 2  The  Merge  Method.  In  1983,  Zicker 
incorporated  the  work  of  Weiss,  Upadhyay,  and  Tenney  (30,31) 
to  form  a  different  method  of  limiting  the  hypothesis 
expansion.  The  purpose  of  this  method  is  to  retain  a  larger 
portion  of  the  entire  measurement  history  (32:26).  Like  the 
"Best  Half  Method",  the  number  of  sequences  is  limited  to  a 
prescribed  depth.  In  this  case,  however,  the  retained 
estimates  are  determined  from  the  probabilistically  weighted 
sum  of  the  estimates  from  hypothesis  sequences  that  differ 
only  in  the  assumption  that  the  oldest  event  in  the  filter 
was  a  signal  or  noise  (See  Figures  2.10  and  2.11).  The 
probabilistic  weight  used  is  the  probability  associated  with 
the  oldest  event  in  the  filter  being  a  signal  or  noise. 
Similarly,  the  retained  variances  are  determined  using  the 
probablistlc  weighting  of  the  oldest  event.  In  equation 


form  this  becomes: 


Hypothesis  sequence 
prior  to  pruning 


Hypothesis  Sequence 
after  pruning 


Snyder -Fishman 
Filter  based 
on  (SS) 


Probability  that  oldest  event 
in  filter  was  a  signal  event 


Snyder-Fishman 
Filter  based 
on  h3  (SN) 


Probability  that  oldest  event 
in  filter  was  a  signal  event 


Snyder-Fishman 
Filter  based 
on  h2  (NS) 


Probability  that  oldest  event 
in  filter  was  a  noise  event 


Snyder-Fishman 
Filter  based 
on  (NN) 


Probability  that  oldest  event 
in  filter  was  a  noise  event 


Figure  2.11  Merge  Method  Pruning  (Depth  *  2) 
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£j,(t)  =  Prtoldest^nJ^CtJ+PrtoldestasJj^tt)  (2-44) 

E^,(t)  »  (Pr (oldest*n)  (P^(t)  +  tx^ft)  -  5Lj,(t)l 

tij(t)  -  2Lj,(t)]T)  +  Pr(oldest=s) (E^ft) 

+  (^(t)  -  Xj,  (t)][^(t)  -  ,  (t )  ]T)  (2-45) 

Note  that  Pr(oldest>n)  +  Pr(oldest=s)  a  i 
where 

2Lj  and  are  the  estimates  associated  with  hypothesis 
that  differ  only  In  the  decision  that  the  oldest 
event  in  the  filter  is  a  noise  or  signal,  respec¬ 
tively.  Here,  j  is  associated  with  the  sequence 
that  assumed  a  noise  event  occurred  and  k  is 
associated  with  the  sequence  that  assumed  a  signal 
event  occurred.  Thus,  k  ■  j  +  2°  1  and 
J  -  0, 1, . . . , ( 2D_1) . 

D  is  the  memory  depth  of  the  hypothesis  tree  before  the 
pruning  process  has  occurred. 

,  and  ,  are  the  "merged"  state  estimates  and  error 
covariances  as  the  number  of  elemental  filters  are 
reduced  from  2°  to  2D-1;  j'  *  0,1, . . ., (20"1) . 

These  are  used  for  propagation  to  the  next 
measurement  time,  hence  the  prime  ( ' )  designation. 

In  this  manner  the  Information  contained  the  both  halves  of 
the  hypothesis  tree  are  retained  to  a  greater  degree  than  in 


the  best  half  method 


One  unique  aspect  of  the  Merge  Method  as  first  noted  by 
Zicker  (32]  is  the  simplification  possible  when  the  Meer 
filter  depth  is  reduced  to  one.  The  Meer  filter 
probabilistic  weighting  Equation  (2-42)  can  be  combined  with 
the  Snyder-Fishman  update  Equation  (2-20)  to  form  a  single 
Meer  filter  update  equation: 

x(t*)  *  x(t“)  +  PrtsJKtt^trj-  H(t1)x(t“)J  (2-46) 

The  Meer  filter  was  designed  to  acheive  the  objective  of 
estimating  the  NPB  centroid  in  the  presence  of  noise.  The 
next  section  will  present  a  MMAE  using  a  bank  of  these  Meer 
filters  to  allow  the  estimation  of  the  NPB  centroid  with  an 
uncertain  beam  time  constant. 

2.6  MMAE  with  Meer  Elemental  Filters 

As  Meer  pointed  out  (21:170),  it  is  important  to  have  an 
accurate  propagation  model  as  the  estimate  of  x(t^)  is  used 
to  approximate  A  by  Equation  (2-7)  with  x(t)  replaced  by 

9 

x(t).  If  the  assumed  model  of  the  underlying  process  has 
characteristics  different  from  those  of  the  actual  physical 
process,  then  the  approximation  of  A_  may  be  Inadequate. 

9 

Since  the  beam  time  constant  is  used  to  propagate  the 
estimate  of  x  via  Equations  (2-15)  and  (2-16)  with 
F  *  — 1/ Tg ,  an  uncertain  beam  time  constant  can  result  in  an 


Inadequate  approximation  of  X  .  Under  these  conditions  we 
we  would  have  degradation  of  performance.  In  other  words, 
the  filter  may  have  poor  robustness.  For  a  filter  variation 
of  only  2%  from  the  true  beam  time  constant,  Johnson  found 
significant  degradation  in  performance  [91.  Johnson 
recommended  HMAE  techniques  be  used  for  online  adaptation  of 
the  beam  time  constant.  The  concept  here  is  to  take  a  set 
of  Heer  filters,  each  modeled  on  a  different  beam  time 
constant,  combine  the  outputs  with  weighting  based  on 
residual  performance  to  arrive  at  a  better  estimate  of  the 
state  and  an  estimate  of  the  actual  beam  time  constant.  The 
MMAE  structure  presented  in  Section  2.4  supports  the 
weighted  estimation  of  the  state.  The  same  structure  also 
supports  the  weighted  estimation  of  the  beam  time  constant. 
If  we  can  estimate  the  correct  beam  time  constant,  then  we 
can  propagate  the  filter  state  estimate  with  greater 
accuracy  by  using  the  beam  time  constant  estimate  in  our 
model.  The  better  state  estimate  provides  a  better  estimate 
of  ^ a  through  Equation  (2-7).  This,  in  turn,  will 
contribute  to  better  state  and  beam  time  constant  estimates. 
K  well  propagated  state  estimate  is  also  important  as  an 
input  to  a  state  feedback  beam  controller.  Similarly,  the 
estimate  of  the  beam  time  constant  could  be  used  in  an 
adaptive  controller  to  determine  the  appropriate  controller 
gains.  Thus,  it  is  desirable  to  seek  an  estimate  of  the 
beam  time  constant.  The  conditional  mean  of  the  beam  time 


constant  Is  given  by  (17:132-133]: 

T(tl)  *  E{^(ti)|Z(ti)=Zi| 

=  e£=1  rk  Pk(t1)  (2-47) 

where  t  is  the  individual  beam  time  constant  associated 
with  each  Meer  filter,  and  pk  is  found  through 
Equation  (2-30).  An  indication  of  the  precision  of  this 
estimate  is  given  by  the  conditional  variance  (17:1331: 

PT  *  E{(  rk-  ^(tj) )2|Z(ti)=Zi  \ 

«  Ek*l(  Tk"  ^(ti))2  Pk(ti)  (2-48) 

The  structure  of  the  MMAE  using  elemental  Meer  filters 
to  estimate  the  beam  state  is  depicted  in  Figure  2.12.  The 
same  structure  used  to  estimate  the  beam  time  constant  is 
depicted  in  Figure  2.13. 

2-7  Summary 

This  chapter  began  with  a  description  of  the  NPB 
dynamics  model.  Following  this,  the  Poisson  space-time 
point  process  description  was  applied  to  the  case  of 
measurements  extracted  from  a  NPB  being  irradiated  with 
laser  beams.  Next,  the  Snyder-Fishman  filter  was  described 
which  uses  the  Poisson  point  process  as  its  measurement 
model.  After  this,  the  MMAE  structure  was  described  to 


support  the  concepts  of  the  Meer  filter  and  MMAE  Meer 
filter.  The  Meer  filter  was  described  as  an  MMAE  structure 
which  uses  a  bank  of  Snyder-Fishman  filters  to  estimate  NPB 
centerline  in  the  presence  of  noise.  Each  Snyder-Fishman 
filter  in  the  bank  operated  under  a  separate  hypothesis  for 
the  order  of  signal  and  noise  events.  Tree  pruning  was  then 
described  as  a  means  to  reduce  the  computational  loading 
caused  by  an  ever-expanding  hypothesis  tree.  Finally,  a 
Multiple  Model  Adaptive  Estimator  (MMAE)  using  elemental 
Meer  filters  was  developed  as  an  approach  to  accommodate  the 
effects  of  an  uncertain  beam  time  constant. 


III.  Controller  Design 


The  MHAE  Meer  filter  provides  the  estimates  necessary  to 
track  the  location  of  the  NPB.  The  next  step  Is  the 
pointing  of  the  NPB  at  a  desired  target.  This  chapter 
develops  a  controller  to  accomplish  that  purpose.  The 
objective  of  the  controller  is  to  point  the  beam  at  the 
target.  Therefore,  we  begin  by  describing  the  target  model. 
Next,  we  will  describe  the  controller  gain  synthesis 
technique.  In  each  of  the  previous  controller  designs, 
assumed  certainty  equivalence  methods  were  employed,  and  the 
controller  gains  were  determined  based  on  the  assumptlo-s  of 
a  linear  system,  deterministic  optimal  control  law, 
full-state  feedback,  and  a  quadratic  cost  criterion.  [8:Sec 
3,1-2;  9:30;  23:22;  32:29].  The  controller  gains  were 
determined  by  LQ  synthesis  techniques  assuming  beam  and 
target  states  were  perfectly  known.  Then,  by  assumed 
certainty  equivalence,  the  full-state  feedback  was  replaced 
by  the  beam  state  estimate  provided  by  the  Meer  filter  and 
the  target  state  estimates  provide  by  a  Kalman  filter. 

Figure  3.1  is  a  simple  block  diagram  of  this  controller 
structure.  This  study  will  use  a  similar  design  process, 
but  will  also  use  Implicit  Model  Following  techniques  to  aid 
the  determination  of  appropriate  controller  gains.  Finally, 
we  will  specify  a  structure  for  the  controller.  In  previous 
thesis  efforts,  the  following  controllers  were  developed: 
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the  proportional  gain  controller,  the  proportional  plus 
integral  (PI)  controller,  and  the  multiple  model  adaptive 
controller  with  elemental  controllers  based  on  separate 
assumptions  of  the  target  dynamics  driving  noise  strength. 
Once  the  controller  structure  is  set,  the  synthesis 
techniques  are  used  to  determine  the  appropriate  gains  to  be 
used  in  the  selected  controller. 
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3 . 1  Target  Models  and  Filter 

We  include  a  target  dynamics  model  in  our  target  state 
estimator  design  to  allow  finer  tracking  of  the  target  than 
possible  with  target  measurements  alone.  This  assumes,  of 
course,  that  the  actual  target  is  well  described  by  the 
target  dynamics  model.  Even  when  the  target  is  well 
described  by  our  dynamics  model,  there  are  uncertalnities  in 
the  dynamics  driving  forces  and  sensor  measurement  data  that 
we  need  to  consider  in  estimating  the  beam  and  target 
locations.  To  combine  the  information  we  have,  or  assume, 
about  the  target  dynamics  and  the  actual  sensor 
measurements,  we  will  use  a  Kalman  filter.  In  incorporating 
a  target  model  into  the  filter,  we  must  consider  the  model's 
impact  on  the  stability  of  the  filter  (and  eventually  on  the 
stability  of  the  closed  loop  system  incorporating  that 
filter).  First,  this  section  describes  the  target  dynamics 
model  used  in  this  study.  Then  the  target  filter  used  to 


3.1.1  Target  Truth  Model.  Our  objective  Is  to  determine 
whether  a  controller,  using  beam  state  estimates  provided  by 
a  Meer  filter,  can  point  the  NPB  at  a  desired  target. 
Therefore  we  need  a  target  model.  In  this  section,  the 
target  truth  model  used  in  this  study  is  presented.  As  this 
is  yet  a  feasibility  study,  we  have  a  wide  latitude  in  the 
choice  of  a  target  model.  Our  desire  is  a  target  which  is 
realistic  and  stresses  the  NPB  controller  performance.  A 
realistic  model  would  include  position,  velocity, 
acceleration  as  functions  of  time,  as  these  are  commonly 
used  to  describe  motion.  As  most  targets  travel  much  more 
slowly  than  the  speed  of  light,  we  will  ignore  relativistic 
effects  in  our  target  model.  A  realistic  model  would  also 
be  three  dimensional,  having  position,  velocity  and 
acceleration  components  in  three  spatial  directions. 

Included  in  the  model  would  be  the  dynamics  driving  forces 
in  each  of  these  directions.  Since  we  have  reduced  the  beam 
model  for  this  study  to  a  single  dimension  in  space,  the 
output  of  a  three  dimensional  target  model  would  have  to  be 
projected  onto  this  one  dimensional  space.  Instead,  we  will 
lmploy  target  model  with  a  single  spatial  dimension  to  match 
the  reduction  of  the  beam  model  to  a  single  spatial 
dimension.  To  account  for  the  order  reduction  we  might 
incorporate  additional  dynamics  driving  noises  on  the  motion 
states  of  the  target  differential  equation  model.  Beginning 
with  Jamerson  in  1985  (8:Sec  3,11],  this  continuing  study 
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has  Incorporated  a  one  dimensional  target  model  having 
components  of  position,  velocity  and  acceleration  with 
dynamics  driving  noise  on  the  acceleration  state.  To 
maintain  continuity  with  efforts  since  that  time,  this  study 
will  use  the  same  model  without  alteration. 

Johnson  describes  this  target  model  as  a  target  dynamics 
position  process  that  is  modelled  as  a  double  Integral  of 
exponentially  time  correlated  acceleration;  a  third-order 
Markov  position  process  that  is  the  output  of  a  third  order, 
linear,  time  invariant  shaping  filter  driven  by  a  stationary 
white  Gaussian  noise  process  that  has  a  mean  of  zero  and 
strength  Q  [9:131.  Figure  3.2  is  a  block  diagram  of  the 
target  model.  The  set  of  first-order  coupled  differential 
equations  relating  the  first-order  Gauss-Markov  acceleration 
to  the  position  are: 


dXTP 

/ 

dt 

— 

> 

X 

(3-1) 

dltTV 

/ 

dt 

= 

XTA 

(3-2) 

dxTA 

/ 

dt 

= 

-l/rT  xta  +  WT 

(3-3) 

or  in  matrix  form: 
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(3-4a) 


( 3-4b) 
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where 


xTp  is  the  target  position  state 
xTV  is  the  target  velocity  state 
x,pA  is  the  target  acceleration  state 

is  the  acceleration  correlation  time  constant  for 
the  target 

wT(t)  is  a  zero  mean  white  Gaussian  noise  of  strength 
QT,i.e.  E jvT(t)vT(t+a)|  =  Qt  fl(a) 

T  subscript  indicates  these  terms  are  related  to  the 
target 

Jamerson  describes  the  position  process  generated  by 
this  model  as  being  more  difficult  to  track  than  a 
first-order  Gauss-Markov  position  process,  but  more 
realistic  [8:Sec  3,11]  for  airborne  targets.  Thus, 
subsequent  thesis  efforts  have  been  based  on  this  target 
model . 

3.1.2  Target  Filter.  Since  one  of  our  objectives  is  to 
determine  the  effects  of  the  Meer  filter  estimates  on  the 
controller,  we  could  measure  controller  performance  with 
target  truth  data  and  Meer  filter  beam  estimates  provided  to 
the  controller.  This  would  be  compared  to  controller 
benchmark  performance  generated  by  providing  the  controller 
with  both  beam  and  target  truth  data.  This  would  isolate 
any  degraded  controller  characteristics  due  to  using  the 
Meer  filter  estimates  associated  with  the  beam.  We  could 
also  provide  the  controller  with  beam  truth  data  and  then 
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estimates  out  of  a  target  Kalman  filter.  This  would  isolate 
any  controller  characteristics  due  to  using  the  target 
filter,  but  development  of  an  optimal  target  filter  is  not 
the  thrust  of  this  study.  However,  there  may  also  be 
degraded  controller  characteristics  that  result  from 
providing  both  beam  and  target  estimates  to  the  controller 
that  do  not  occur  with  beam  estimate  and  target  truth  data. 
So  that  we  can  determine  controller  characteristics  when 
Meer  filter  estimates  are  combined  with  target  filter 
estimates,  i.e.,  so  that  we  can  evaluate  the  performance 
potential  of  the  eventual  implementable  algorithm,  we  must 
design  a  target  filter.  This  section  presents  the  target 
filter  that  provides  the  target  estimates  to  the  controller. 

In  designing  our  target  filter  we  would  like  to  exploit 
all  the  information  we  have  concerning  the  target,  our 
assumptions  about  the  driving  dynamics,  and  our  assumptions 
about  measurement  errors  committed  by  our  sensor  system.  To 
combine  this  Information  optimally,  we  will  use  a  Kalman 
filter.  As  the  focus  of  this  study  is  on  the  Meer  filter, 
it  is  reasonable  that  the  target  filter  model  should  be 
based  on  our  truth  model.  The  forms  of  the  dynamics  model 
and  the  measurement  model  for  the  target  are  linear 
time-invariant  stochastic  equations  driven  by  white  Gaussian 
noise,  so  the  standard  Kalman  filter  equations  can  be  used 
to  determine  the  target  state  estimates.  Therefore,  the 
Kalman  filter  state  estimate  and  error  covariance  matrix 
propagation  equations  become: 
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2Sup ^  ^  i + 1 )  = 


VW  *  4T(ti+l,ti)ET(ti,4T<ti+l,ti 


♦T(tii.i^«)G(<*)Q(«)GT(«)*^(ti+i/<*)  <3« 


(3-5) 


(3-6) 


where  is  the  state  transition  matrix  associated  with  F,f: 


=  I  o 


At  T*lAt/TT  -  1  +  exp(-4t/rT)] 
1  rT ( 1-exp ( -4t/rT ) 1 

0  exp(-4t/rT) 


(3-7) 


where  At  Is  the  sample  period. 

Ve  have  further  assumed  that  measurements  of  the  target 
position  can  be  made  available  at  regular  intervals  rather 
than  the  random  intervals  as  with  the  Meer  filter.  For  this 
study,  a  sensor  such  as  a  radar  or  electro-optical  sensor  is 
assumed  to  provide  a  direct  measurement  of  the  target 
position  states,  with  some  measurement  corruption  noise 
added  to  account  for  nonlinear  effects  and  variations  in  the 
index  of  refraction  in  the  atmosphere  (for  near  earth 
targets),  sensor  vibration,  background  radiation, 
non-signal-induced  currents  in  the  sensor,  and  fundamental 
limits  due  to  the  operating  wavelength  of  the  sensor.  The 
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Central  Limit  Theorem  supports  the  fact  that  the  combined 
effects  of  these  noise  sources  may  be  modelled  as  Gaussian. 
Therefore,  we  will  model  the  noise  as  a  zero  mean  white 
Gaussian  noise  with  variance  R.  For  simplicity,  it  is  also 
assumed  that  these  measurements  can  be  expressed 
conveniently  in  the  same  coordinate  frame  as  the  beam 
location.  To  account  for  this  explicitly  would  merely 
involve  a  transformation  matrix.  This  results  in  a 
measurement  equation  of  the  form: 

z(ti)  =  H^tj)  xT(t1)  +  v(t^)  ( 3-8a ) 

or  in  this  case, 

zttj)  =[1001  xT(tl)  +  v(t1)  ( 3-8b) 

where  z  is  a  scalar  measurement  and  v  is  assumed  to  be  a 
zero-mean  white  Gaussian  discrete  time  noise  with  variance 
R.  The  measurement  noise  v  is  assumed  to  be  independent  of 
the  dynamics  driving  noise  w  of  the  target.  H  reflects  our 
assumption  that  the  sensor  measures  the  target  position 
state  directly.  Figure  3.3  depicts  the  target  and  sensor 
models  upon  which  the  filter  is  based. 

With  these  assumptions  the  standard  Kalman  filter 
measurement  update  process  is  given  by: 

iyCtJ)  -  VV  +  jy^Mztt^  -  Uptt^iytt ~)1  (3-9) 
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EjUj)  *=  Er(tj)  -  K{ti)HT(t1)ET(t^)] 


(3-10) 


K(tL)  =  PT(t")iLj(ti)  (H^t^p^t")  h£<V  +  a(  t± ) ) -1 


(3-11) 

A  A 

with  initial  conditions  x(tQ)=  Xq  and  P(tQ)  =  Pg. 

Figure  3.4  is  the  Kalman  filter  block  diagram.  Note  that 
there  is  no  feedback  to  the  target,  and  so  there  is  no 
feedback  control  term  in  Equations  (3-5)  or  (3-9),  or  in  the 
filter  diagram.  The  state  estimates  out  of  this  filter 
provide  the  target  state  estimates  used  by  the  controller. 

Before  we  incorporate  the  target  model  and  sensor  model 
into  the  Kalman  filter,  we  should  verify  that  the  resulting 
filter  is  stable.  Therefore,  we  must  examine  the  filter 
model's  stabllizability  with  respect  to  the  dynamics  driving 
noise  input  w  and  detectability  with  respect  to  the 
measurement  z  output  of  the  sensor  system  (12:464-465; 
16:243-2451.  By  stabllizability,  we  mean  the  filter's 
unstable  subspace  is  contained  in  its  controllable  subspace 
(12:4621.  By  detectability  we  mean,  the  filter's 
unreconstructible  subspace  is  contained  in  its  stable 
subspace  (12:4651.  As  noted  by  Maybeck  (18:92-931:  "By  its 
definition,  stabllizability  is  implied  by  either  complete 
controllability  or  asymptotic  stability...  It  can  be 
readily  be  seen  that  either  complete  observability  or 
asymptotic  stability  will  imply  detectability."  Since  the 
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conditions  for  controllability  and  observability  are 


stricter  than  stabi lizabillty  and  detectability, 
establishing  controllability  and  observability  will 
establish  the  stability  of  the  Kalman  filter.  For  linear 
systems,  the  conditions  for  asymptotic  stability  are 
indistinguishable  from  conditions  for  bounded  input-bounded 
output  (BIBO)  stability.  Thus,  showing  that  either  set  of 
conditions  is  satisfied  will  indicate  that  the  filter  is 
asymptotically,  globally  stable  (16). 

The  system  model  described  in  Equation  (3-4)  is  said  to 
be  stochastically  controllable  if  there  exists  positive 
numbers  a  and  ft,  0  <  a  <  ft  <  «,  and  a  positive  integer  N 
such  that,  for  all  i£N, 


«I  < 


j=i-N+l 


*(ti,tj)Gd(tj_1)Qd(tj_1)gJ(tj_1)*T(ti,tj)  <  ftl 

(3-12) 


The  system  model  is  said  to  be  stochastically  observable 
if  there  exists  positive  numbers  a  and  ft,  0  <  a  <  ft  <  and 
a  positive  integer  N  such  that,  for  all  i£N, 


al  5  £  ♦T(tj,t1)HT(tj)g.‘1(tj)H(tj)i(tJ,t1)  £  ftl  (3-13) 

3=i-N+l 

The  chosen  models  do  result  in  observable  and 
controllable  system  models.  Therefore,  for  our  choice  of 
dynamics  and  sensor  models,  the  target  filter  is  guaranteed 
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to  be  stable.  In  a  similar  manner,  we  will  next  chose  and 
verify  the  model  to  be  used  in  the  controller  development. 

3 . 2  Controller  Gain  Via  LQ  Synthesis. 

With  the  controller  structure  specified  as  in  Figure  3.1 
(more  will  be  said  about  this  in  Section  3.3),  we  need  a 
method  of  determining  the  appropriate  feedback  gains.  In 
this  study  we  will  use  LQ  synthesis  techniques  to  generate 
the  gain  appropriate  to  a  deterministic  full-state  feedback 
controller.  Then,  by  applying  assumed  certainity 
equivalence  we  will  use  these  gains  for  the  NPB  controller. 
Describing  LQ  synthesis  will  involve  explaining  the 
quadratic  cost  function.  The  solution  to  the  minimization 
of  the  quadratic  cost  function  results  in  a  method  of 
determining  the  optimal  gain.  Since  the  cost  function  and 
therefor**  the  gain  involve  choosing  weighting  matrices,  we 
will  present  Implicit  Model  following  as  an  aid  to 
determining  appropriate  weighting  matrices.  This  section 
will  discuss  the  use  of  the  quadratic  cost  to  aid  the 
determination  of  controller  gains. 

3.2.1  Description  of  Quadratic  Cost  Function.  The 
purpose  of  the  cost  function  is  to  represent  what  we  want  to 
happen  to  the  system  when  it  deviates  from  our  desired 
response,  or  how  we  want  the  system  to  be  controlled,  so  an 
approach  to  controller  design  can  be  developed.  To  see  how 
the  system  model  is  incorporated  into  the  weighting  matrices 
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of  the  cost  function,  consider  the  system  described  by  the 
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stochastic  differential  equation. 


2L  ( t )  =  F  ( t )  x  ( t )  +  B(t  )u(  t )  +  G(t)w(t) 


e{  da(t)  |  =  0 

(  T  )  t-f 

E<  dfl(t)  da  (t* )  J  =  ; 


(3-14) 


or  more  precisely, 

dx(t )  =  F(t)x(t)dt  +  B(t)u(t)dt  +  G(t)da(t)  (3-15) 

where : 

a  is  a  Brownian  motion  process  with  statistics 


( 3-16a ) 

( 3-16b ) 


& 


0  otherwise 


£ 


« 


The  solution  to  Equations  (3-14)  or  (3-15)  is: 


r." 


>s 


f* 


x(t)  *  ♦(t,t())x(t0)  £(t,a)B(a)u 


(a  )da 


P. 

J  40 


(t,a)  G(a)dB(flt) 


(3-17) 


and  the  equivalent  discrete-time  system  is  described  by  the 
stochastic  difference  equation: 

-(ti+l)  =  ft(ti+l,ti)-(ti)  +  &c|(  t±  )u< i )  +  G^t^WjU^ 

(3-18) 


where : 


£  is  the  state  transition  matrix  associated  with  £ 

is  the  discrete  time  input  distribution  matrix  given 
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(3-19) 


and  u(a)  Is  assumed  to  equal  u(t^)  for  all  a  In 
the  interval  lti,ti+l) 

Gj  is  the  discrete  time  driving  noise  distribution 
matrix, 

equal  to  the  identity  matrix  if  Q.  is  evaluated  as 
done  below 


Wj  is  an  equivalent  white  noise  with  statistics: 


E  1  Brt  1  -  0 


e!  )  = 


( 3-20a ) 

°d(ti)  fci=tj 

( 3-20b) 

0  otherwise 


where 


i+1 


t(ti+i/«)fi(«)ft(«)ST(«)iT<ti  +  1/«)  <J«  (3-21) 


Using  the  LQG  assumptions  and  applying  assumed  certainty 
equivalence,  we  will  use  LQ  synthesis  techniques  to 
determine  the  controller  gains.  Therefore  we  can  drop  the 
stochastic  term  in  the  subsequent  development.  One  approach 
to  representing  the  system  in  a  cost  function  is  applying 
quadratic  penalties  to  deviations  in  the  states  and  controls 
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from  zero.  Then  the  cost  function  has  the  following  form 
for  a  sample  data  controller  (18:731: 


J  «=  E(  l/2jlT(tN+1)XF2L(tN+1)  +  E  l/2[2«T(t1)2L(t1)2.(t1) 

♦  aT(t1)U(t1)Ji(t1)  +  22cT(t1)£(t1)y.(t1)l  )  ( 3-22a ) 


J  -  E(  l/2xT(tN  +  1>3CFx(t|I+1) 


+  Z  1/2 
1=0 


)~|T  rac<ti )  s.(t1nr2c(t1) 

.uttjlJ  L  £T  <  1 1 )  U<  t  L  )  JLsi<  > 


( 3-22b ) 


where ; 


E  is  the  expectation  operator 
tQ  is  the  initial  time 

fcN+l  ls  the  flnal  time 

tN  is  the  last  time  control  is  computed  and  applied  as 
a  constant  over  the  ensuing  sample  period 
N  is  the  number  of  sample  periods  from  tQ  to  the  final 
time  of  control  application  t^ 
x  ls  the  controller  state  vector, 
u  is  the  control  input  vector. 

X_,  Hr  \Lr  and  S  are  weighting  matrices  described  below. 


Xp  and  X  are  symmetric  positive  semldefinite  weighting 
matrices  used  to  emphasize  the  importance  of  minimizing  each 
system  state  deviation  from  zero  (perturbation  states  are 
driven  towards  zero)  at  the  final  time  and  each  incremental 
time,  respectively.  A  large  value  for  a  diagonal  element  of 
JC,  as  compared  to  other  diagonal  elements,  would  result  in 
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the  design  of  controller  gains  that  would  maintain  the 
system  state  associated  with  that  diagonal  element  closer  to 
zero  than  the  other  system  states. 


U  is  a  symmetric  positive  definite  matrix  used  to  weight 
the  importance  of  minimizing  control  application.  A  smaller 
value  for  a  diagonal  element  of  U,  as  compared  to  other 
diagonal  elements,  would  result  in  controller  gains  that 
would  allow  more  energy  to  be  used  in  the  actuator 
associated  with  this  diagonal  element  than  allowed  through 
other  actuators.  The  choice  of  elements  in  X  and  the 
corresponding  entries  in  U  forms  a  tradeoff  between  the 
desire  for  minimum  rms  system  states  and  minimum  control 
energy  application. 

S  is  a  matrix  that  arises  naturally  from  a 
continuous-time  cost  function.  Given  a  system  described  by: 

x(t)  =  F(t)x(t)  +  B(t)u(t)  (3-23) 


with  a  continuous  time  quadratic  cost  of: 

V  E  1/a*T' WM'W 


X(t  )' 

T 

Xx(t)  Sxu(t>' 

JL(t)' 

dt 

fto 

L«xu  «uu(t,J 

,U(t). 

(3-24) 


where  is  a  positive  semidefinlte  state  weighting  matrix, 

S^u  is  a  positive  definite  control  input  weighting  matrix 

and  w  is  a  cross  weighting  matrix.  Note  that  these  are 
— xu 

continuous  time  weighting  matrices.  It  can  be  shown 
(18:73-751  that  the  weighting  matrices  for  the  equivalent 


discrete  time  system  are  given  by/ 


.•  •/:: 


t  i  (t/t1)S?xj{(t)*(t,t1)]dt 


( 3-25) 


mft. 


i  i  (t/t  >w  (tjitt/t.)  +  w  (t) 

i  X  -UU  (3-26) 

+  IT(t/t1)Wxu(t)  +  S^u(t)|(t,ti)  Idt 


Tti+1  T 

Sftj)  »  I  (  i^t/tj 
t  a 


JWj^ftJlft/tj)  ♦  ♦(t/t1)wxu(t) )dt 

(3-27) 


where 


:/tt)  *  *(t, «)&(<*; 


(3-28) 


Thus,  the  S  terms  arises  from  two  sources.  The  first  is  a 
desire  for  cross  terms  in  the  discrete  time  representation, 
i.  e.  nonzero  W^lti  terms  above,  or  for  placing  quadratic 
penalty  on  deviations  of  x(t)  from  zero  (See  Equation 
(3-23)).  The  second  is  a  desire  to  exert  control  influence 
on  the  states  over  the  entire  sample  period,  i.  e.,  the  part 
of  Equation  (3-27)  other  than  that  due  to  nonzero  tf^ft). 
Also  S  is  chosen  so: 


i(tl)  s.<v 

^(tj)  LL( tj ) 


is  positive  semidefinite  for  all  t*.  This  is  done  to  ensure 


that  the  solution  to  the  discrete  backward  Riccati  recursion 
converges  to  a  constant  steady  state  solution  [18:95]. 

Use  o£  a  quadratic  cost  for  the  control  state  terms  (x) 
in  Equations  (3-228)/  ( 3 — 2 2b )  and  (3-24)  arises  from  the 
desire  to  minimize  the  mean  square  error.  The  square  error 
is  useful  as  a  cost  basis  as  it  provides  for  very  large 
increases  in  cost  as  deviations  in  the  control  state  from 
zero  grow,  while  allowing  for  smaller  increases  in  cost  as 
deviations  in  the  control  state  remain  near  zero.  Use  of 
the  square  of  the  control  states  also  prevents  any 
deviations  to  result  in  a  lower  or  negative  cost. 

Quadratics  are  also  appropriate  for  the  control  energy. 

Use  of  a  quadratic  cost  for  the  control  energy  terms  (u) 
in  Equations  (3-22a),  (3-22b)/  and  (3-24)  arises  naturally 
from  the  common  form  of  energy  equations.  Energy  equations 
come  in  the  form  1/2  mv‘  for  a  mass  with  translational 

velocity  v  or  1/2  Kx2  for  a  spring  with  displacement  x  or 

9  2 

1/2  Cv*  for  a  capacitor  with  voltage  v  or  1/2L1  for  an 

inductor  with  current  i  [3: 31; 181.  Thus,  a  general  control 

energy  can  be  expressed  in  terms  of  a  coefficient  times  the 

control  variable  squared. 

As  Naybeck  describes  the  use  of  this  cost  function  with 

linear  systems  and  Gaussian  noise  models  [16:12]: 

Although  other  types  of  cost  functions  are  useful, 
quadratics  are  a  good  physical  description  of  many 
control  objectives,  such  as  minimizing  mean  square 
errors  while  not  expending  inordinate  amounts  of 
control  energy.  Moreover,  if  one  views  most 
linear  systems  as  being  derived  as  linear 
perturbation  models,  by  developing  a  controller 


that  specifically  attempts  to  maintain  small 
values  of  quadratics  in  state  and  control 
deviations/  one  is  inherently  enhancing  the 
adequacy  of  the  linear  perturbation  model  itself. 
However /  perhaps  the  most  important  reason  is  that 
this  combination  of  modelling  assumptions  yields  a 
tractable  problem  whose  solution  is  in  the  form  of 

y  implemented 


feedback  control  law. 


Once  the  quadratic  cost  function  is  established/  we  will 
seek  the  solution  that  results  in  the  minimum  cost.  The 
next  section  will  present  the  solution  and  show  how  the 
optimal  gain  is  determined  from  this  solution. 


3.2.2  Controller  Gain  Synthesis.  Equipped  with  Linear 
models  and  a  Quadratic  cost  function  we  can  synthesize  an 
optimal  linear  full-state  feedback  control  law.  As  an 
example,  this  will  be  done  for  a  Proportional  plus  Integral 
(PI)  controller.  Details  regarding  the  PI  controller  are  in 
Sections  3.3.3  and  3.3.4.  The  gain  for  a  proportional  gain 
controller  (regulator)  is  determined  in  the  same  manner. 

First  we  consider  an  augmented  system  (see  Section  3.3.3 
for  further  discussion  of  the  augmented  system)  described  by 
the  state  equation  which  is  both  controllable  and 
observable : 


^a(ti+l)  =  ®a(ti+l'tl)  W  +  (3-29) 


The  optimal  control  law  is  118:69): 


a  (tj)  =  «■  ntj)  % 


<3-30a) 
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where 
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t 

* 

gc(ti>  is  the  optimal  controller  gain 

q  is  a  pseudointegral  state  described  in  Sections  3.3.3 
and  3.3.4 

yd  is  the  nonzero  setpoint 

E(t^)  is  the  gain  applied  to  achieve  the  desired 
nonzero  setpoint 


Defining  optimality  as  minimizing  the  quadratic  cost 
function  presented  in  the  last  section,  then  that  minimum 


cost  is  given  by: 


l/2[jcA(t1)Kc(ti)x(ti)  ) 


(3-31) 


where  K  (t.)  is  the  solution  to  the  backward  Riccati 
— c  1 

T  T 

equation  for  the  augmented  state  vector  [x  q]  : 


MV  *  l*WV  -  MV  MV1  MVi* 
‘♦<Vi'V  "  MV  MV1  +  MV  U(V  MV 


♦  i(ti)  -  attj)  ^(t^  -  g*:(tL)’T  a(ti)T 


(3-32) 


MW  * 


(3-33) 


where  £  is  the  state  transition  matrix  of  the  controller- 
assumed  augmented  model  for  gain  calculations  and  is  the 
corresponding  control  distribution  matrix.  The  optimal 


iUhSUV 


controller  gain  Is  given  by: 


S^V  -  fii(ti)+iTj(ti)!(c(ti+1)Bd(t i)T1 

l8d(tl,Sc(tiil)i(tl+l'ti>+iT(tl)1  (3-34) 

We  will  seek  to  develop  a  constant  gain  controller,  as 
this  will  be  the  least  computationally  burdensome  and  it 
will  not  impose  unacceptable  performance  degradation.  To 
accomplish  this,  the  weighting  matrices  X,  U,  and  S  will  be 
assumed  to  be  time  invariant,  and  the  state  transition 
matrix  $.(ti4l,ti)  and  control  distribution  matrix  gg  are 
also  time  invariant.  In  addition,  by  letting  N  approach 
infinity  in  Equations  (3-22a)  and  3-22b),  the  effect  of  the 
cost  at  the  time  t  N+i  on  the  total  cost  is  assumed  to  be 
negligible.  In  other  words,  the  terminal  transient  in  the 
Riccati  equation  is  short  compared  the  the  time  for  the 
entire  process  for  this  application.  This  reduces  the 


equations  to: 


(3-35) 


(3-36) 


u  (t^  =  -g^tj) 

s£  -  <im5W1  <bX.i*4t> 


where  is  the  solution  to  the  steady-state  backward 
Riccati  equation  118:72-76], 


Kc  -  (3-37) 


Similarly,  a  constant  gain  E  for  Equation  (3-30)  is  found  to 
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be : 


E  =  (Scl  -  Gc2  E^12/Kc22  l!L12  +  *22 


where 
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(3-39) 
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(3-40) 


Thus  we  have  a  method  of  determining  the  controller 
design  based  on  the  weighting  matrices  we  choose.  In  the 
next  section,  the  use  of  Implicit  Model  Following  is  used  to 
determine  a  set  of  weighting  matrices  for  use  in  the 
quadratic  cost  function. 

3.2.3  Weighting  Matrices  Using  Implicit  Model  Following. 
Model  following  in  general  is  a  means  of  forcing  the  outputs 
of  a  controlled  system  to  behave  as,  or  mimic,  a  desired 
system  model  that  meets  specifications  [22:Sec  11,2).  There 
are  two  types  of  model  following  controllers:  explicit  model 
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followers  and  implicit  model  followers  [22:Sec  II, 2].  This 
section  will  briefly  describe  the  explicit  model  follower 
(as  it  is  not  used  in  this  thesis)  and  then  describe  the 
implicit  model  follower.  It  is  the  implicit  model  follower 
that  is  used  as  an  aid  to  determining  the  weighting  matrices 
in  this  effort.  These  weighting  matrices  form  the  cost 
function  and  are  used  to  determine  controller  gains  via  the 
backward  Riccati  equation  (Equations  (3-32),  (3-33),  and 
(3-34),  or  (3-36)  and  (3-37)  for  the  steady  state  constant 
gain  controller). 

3.2. 3.1  Explicit  Model  Following.  In  explicit 
model  following,  the  desired  model  is  explicitly  simulated 
in  the  controller.  The  outputs  of  the  desired  model  are 
multiplied  by  feedforward  gains  that  have  been  generated  by 
the  controller.  These  control  inputs  drive  the  system  to 
respond  as  the  explicit  model.  Thus,  the  explicit  model 
follower  requires  additional  online  computer  processing  to 
provide  the  desired  model  outputs.  In  general,  the  desired 
model  is  composed  of  simple  first  or  second  order  systems 
for  each  response  mode  of  interest,  so  the  computational 
burden  may  be  acceptable.  Except  that  the  target  might  be 
loosely  considered  an  explicit  model,  this  study  does  not 
incorporate  an  explicit  model.  Therefore,  explicit  model 
following  will  not  be  developed  further. 

3. 2. 3. 2  Implicit  Model  Following.  In  Implicit 


Model  Following,  the  desired  model  dynamics  are  introduced 


into  the  controller  gain  synthesis  through  the  cost  function 
via  the  weighting  matrices.  There  are  no  model  states  to 
follow,  so  model  following  is  only  implied,  hence  the  name 
implicit  model  following.  The  controller  designed  through 
implicit  model  following  provides  feedback  whenever  the 
system  dynamic  response  deviates  from  the  desired  model 
dynamic  response.  As  such,  implicit  model  following 
attempts  to  match  the  system's  closed-loop  poles  with  those 
of  the  desired  model.  This  ties  the  closed-loop  system 
characteristics,  including  stability  robustness,  to  the 
desired  model  characteristics.  (18:Sec  11,8).  As  Miller 
points  out  (22:Sec  V,15),  there  is  a  limit  to  the  number  of 
eigenvalues  and  eigenvectors  that  can  be  independently 
assigned  in  the  closed-loop  system  by  implicit  model 
following.  This  limit  is  defined  by  the  number  of 
independent  outputs  and  controls  of  the  system.  In  Miller's 
words  (22: Sec  V, 151: 

If  the  controlled  system  has  r  independent 
controls  and  p  Independent  outputs  available  for 
feedback,  then  at  most  max(r,p)  eigenvalues  may  be 
assigned  and  min(r,p)  entries  of  max(r,p) 
closed-loop  eigenvectors  may  be  assigned. 

For  our  NPB  controller,  we  have  only  one  independent  control 
and  one  independent  output.  Thus,  only  one  eigenvalue  and 
one  eigenvector  may  be  Independently  assigned  in  the 
closed-loop  system.  Through  Implicit  Model  Following  we  can 
drive  the  system  toward  the  desired  pole.  In  Miller's 


words,  referencing  Krelndler  and  Rothschild  (22:Sec  11,6-7; 
11:835-8421 : 

The  implicit  formulation  has  an  effect  that  is 
asymptotically  (as  state  weighting  grows  large) 
equivalent  to  eigenstructure  assignment,  and  the 
matching  of  the  poles  of  the  model  and  controlled 
system  can  produce  desirable  transient  response. 
Rejection  of  unmodelled  zero-mean  disturbances  may 
also  be  better  than  with  an  explicit  controller, 
since  good  disturbance  rejection  characteristics 
can  be  included  in  the  model  response  so  that  the 
disturbance  states  need  not  be  modelled  in  the 
controller . 

There  are  a  number  of  ways  to  incorporate  the  desired 
model  into  the  cost  function.  Maybeck,  Miller,  and  Howey 
[151  presented  a  method  of  defining  the  discrete  time 
weighting  matrices  from  the  continuous  time  dynamics  model, 
continuous  time  weighting  matrices  and  continuous  time 
desired  dynamics  model.  Johnson  (9:173-1801  developed  an 
alternate  approach  that  includes  the  state  transition 
matrix.  Since  we  are  developing  a  discrete  time  controller 
using  discrete  time  system  models,  we  may  desire  to  begin 
with  discrete  time  weighting  matrices.  Therefore,  we  will 
present  a  method  for  determining  the  appropriate  implicit 
model  following  weighting  matrices  from  the  discrete  time 
system  description  and  a  set  of  discrete  time  weighting 
matrices . 

To  see  how  the  desired  model  is  incorporated  into  the 
weighting  matrices  of  the  cost  function,  consider  the  system 
described  by  Equations  (3-14)  to  (3-21).  Again,  we  will  be 
using  LQG  assumptions  and  applying  assumed  certainty 
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equivalence  so  that  we  can  use  LQ  synthesis  techniques  to 


determine  the  controller  gains.  As  before,  we  will  drop  the 


stochastic  portion  of  the  models  to  be  used.  The  sample 


data  equilvalent  cost  function  without  the  incorporation  of 


the  implicit  model  is  given  by  Equation  (3-22a): 


J  -  s{  l/2l!t(t(t+1)'VtN11)(![(tN(1)) 


n  _  _ 

+  (  1/2 (  xT(ti)X(t1)x(ti)  +  uT(t1)LI(ti)u(ti) 


+  2xT<t1)g.(t1)ii(t1)  )  ]}  ( 3-22a ) 


Our  real  interest  is  in  the  continuous  time  outputs  of  the 


system  is  given  by. 


y.(t)  =  C(t)£(t) 


(3-41) 


where  C(t)  relates  the  controlled  outputs  to  the  system 


states.  However,  to  the  sample  data  controller  the  outputs 


of  interest  are. 


*  £(t1)jt(t1) 


(3-42) 


The  objective  is  to  force  the  output  dynamics  to  match,  as 


closely  as  possible,  the  output  dynamics  of  the  desired 


model,  described  by  the  equations: 


^(t)  «  ^(t) 


(3-43) 


or  in  the  discrete  time  representation 


^m^i+l*  *  4m(ti+l'ti)*m(ti) 


(3-44) 
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where  ym  is  of  the  same  dimension  as  y  and  £m  is  the  state 
transition  matrix  associated  with  the  desired  dynamics  model 
Fm.  This  model  describes  a  closed-loop  system  which  is 
primarily  robust,  but  also  has  good  performance 
characteristics.  An  appropriate  £ncan  be  determined  from 
the  desired  eigenstructure,  with  poles  placed  in  desired 
locations  and  eigenvectors  as  nearly  orthogonal  as  possible 
for  robustness  [71. 

To  design  a  controller  using  a  continuous  cost  function 
and  implicit  model-following,  we  would  place  a  quadratic 
penalty  on  deviations  of  the  continuous  time  output  from  the 
desired  characteristics  or  (y(t)  -  F  (t)y(t)J.  Analogously, 
we  will  place  a  quadratic  penalty  on 
yft^ll  with  a  symmetric  positive  semidefinite  output 
weighting  matrix  X j  ( t A ) .  The  cost  function  now  takes  the 
form: 

J  .  e{  VW 

N  T 

♦  E  l  l/2(ty(ti+1)-tm(t1+1,t1)y(ti)J1  XJ(tl) 

Iz(ti+l)‘4m(ti+l'ti )x(ti) 1 

♦  2l^(ti+i)-ftm(ti+i#ti)z(ti> lT(ti)s-(ti,tt(ti)  J  }  (3-45) 

By  applying  Equations  (3-23)  and  (3-42),  the  cost  function 
can  then  be  rewritten  as: 


‘ | d IS' J*»  I* k.1  #,%  h1Ji,Ji*A',A*.H*»i fc  .4  . 


N 

+  E  (  1/2 (  2tr(t1)Jt1(t1)Jt(t1)  +  Jl  (t1)UI(t1)ji(t1) 


♦  2x 


xT(ti)£I(t1)ii(t1)  ]  |  (3-< 


where  (dropping  time  index  for  clarity): 


*1  - 

&!  -  IS*-tn£]TII£fid  +  (£♦-*„£)  TS. 


(3-47) 


(3-48) 


Uj  =  (£^jlTXI£Bd  +  C  1TS.  +  STCBj+  U  (3-49) 


where  the  I  supscript  indicates  an  implicit  model  following 
weighting  matrix.  Taking  the  cost  function  summation  out  to 
infinity  for  the  steady  state  condition. 


I1 

•  pLtti)]*  [VV  £i(ti)l  fJt(ti)l 

1/2  (  E  )  (3-50) 

i-o  ll!(ti>J  La(tl)J 

Equations  (3-46)  and  (3-50)  are  in  the  familiar  cost 
function  form,  for  which  we  can  determine  the  controller 
gains,  using  the  backward  Riccatl  equation  as  in  Section 
3.2. 

The  form  of  the  cost  function  shown  allows  the  use  of 
standard  LQ  synthesis  techniques  to  develop  an  optimal 
feedback  controller  that  tries  to  force  the  system  to  meet 
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classical  performance  specifications  embodied  in  the  model 
output,  y  .  The  extent  to  which  the  desired  response  is 
achievable  is  dependent  on  the  degree  of  difference  between 
the  system's  inherent  dynamics,  represented  by  and  the 
desired  response  dynamics,  represented  by  *m  (22:Sec  11,3-5] 
and  how  heavily  ly(ti+1)  “  ♦  vt t ^)  1  is  weighted  through 
Equations  (3-47),  and  (3-48).  The  desired  model,  repre¬ 
sented  by  ♦m,  is  the  goal,  but  it  is  the  weighting  on  Xj 
that  determines  how  much  control  energy  is  applied  to 
achieve  this  goal.  This  form  also  allows  the  experience 
gained  in  an  iterative  design  process  to  be  useful  when 
applying  Implicit  Model  Following  techniques. 

The  choice  of  an  implicit  model  must  be  made  carefully. 
Since  the  feedback  gain  is  a  function  of  the  difference  of 
the  system  model  and  the  desired  model  through  the  matrix 
(Ot-^C),  the  Implicit  model-follower  is  part icular ily 
sensitive  to  parameter  variations  and  plant  modelling 
inaccuracies  (22:Sec  11,71.  In  this  study  we  are  already 
aware  of  a  problem  with  the  mismodelling  of  the  beam  time 
constant.  Therefore,  we  should  be  concerned  about  using 
implicit  model  following.  However,  if  elements  in  *m 
dominate  elements  in  £,  then  for  small  variations  in  £., 

(£*-#  C)  will  remain  essentially  the  same.  Even  so,  the 
elements  of  should  not  be  set  so  as  to  result  in  an 
unstable  clcsed-loop  system.  As  such,  there  may  be 
limitations  to  the  performance  or  robustness  Improvement 
possible  when  using  implicit  model  following. 


As  with  the  selection  o£  the  target  model,  we  have  a 
great  deal  o£  freedom  in  chosing  a  controller.  In  chosing  a 
controller,  we  desire  minimum  transient  times  and  minimum 
steady  state  errors  £or  wide  variations  in  the  parameters  of 
the  closed-loop  system.  Also,  we  desire  a  computationally 
fast  structure.  These  desires  shape  the  complexity  of  the 
controller.  This  section  presents  the  controller  used  in 
this  thesis.  Prior  to  discussing  the  controller  structure, 
some  of  the  basic  controller  building  blocks:  Proportional 
Gain  (PG)  regulator,  PG  tracker,  PI  regulator,  PI  tracker, 
and  Multiple  Model  Adaptive  Controller  (MMAC)  are  presented. 
Then  we  will  present  the  controller  used  in  this  study. 

3.3.1  Proportional  Gain  Regulator.  The  purpose  of  the 
regulator  is  to  drive  the  beam  centroid  to  the  center  of  the 
detector  array.  The  discrete-time  state  equation  of  the 
controlled  beam  is  given  by: 


VW  *  VWWV  *  Wulti>  ♦  ww 

(3-51) 

The  optimal  control  law  is  given  by: 


u*(ti)  =  -G*{ti JXgttj) 


(3-52) 


The  optimal  control  gain  is  calculated  as  is  Section  3.2.2 
with  the  gain  E  being  ignored.  Since  we  are  interested  in 
directing  the  beam  onto  the  target,  we  will  consider  the 


tracker  next. 


3.3.2  Proportional  Gain  Tracker.  Maybeck  and  others 
have  demonstrated  how  LQG  techniques  can  be  used  to  develop 
trackers  (6:114-1221.  LQ  synthesis  techniques  require  that 
the  model  used  to  determine  controller  gains  have  controll¬ 
ability  with  respect  to  the  control  input  u  and  observ¬ 
ability  with  respect  to  the  commanded  output  y  .  For  the 

c 

time-invariant  case,  the  controller  needs  stabilizability 
with  respect  to  the  control  feedback  u  and  detectability 
with  respect  to  extracting  the  controlled  output  y  .  This 

v 

is  sufficient  to  ensure  a  solution  to  the  backward  Riccati 
equation.  Using  LQG  techniques  and  assumed  certainty 
equivalence,  we  first  consider  the  deterministlstic  LQ 
optimal  full-state  feedback  controller  based  on  the  beam  and 
target  models  with  the  uncertainities  removed. 


2L<t  )  -  1  (t  ,t  )£(t  )  +  £  (t  )u(t  )  (3-53) 

i+1  C  i+1  1  i  dC  i  i 


where 


♦b  I  0 


2-  1  *T- 


(3-54) 


and  & 


IB  0  0  0 IT,  with  B  given  by  Equation  (3-19). 


The  controller  continuous-time  dynamics  model  is  described 


in  terms  of: 


PB  I  0 

1  |  E. 


(3-55) 
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For  the  tracker  we  desire  to  point  the  beam  at  the  target. 
Another  way  of  expressing  this  desire  is  to  drive  the 
difference  between  the  beam  output  and  the  target  postion  to 
zero.  Then  deviations  from  zero  can  be  considered  as  errors 
given  by: 

e(t)  =  xB(t)  -  xTp(t)  (3-56) 

where  the  target  position  is  taken  as  the  reference 
variable.  These  errors  also  define  the  output  state  of  the 
controller . 

yc  *  cc£  -  tcB  Ofix 

yc  «  U  -1  0  °J*  (3-57) 

therefore  y  is  the  variable  we  want  regulated  to  zero. 

Again,  we  need  controlability  or  at  least 
stabilizability  with  respect  to  the  control  input  u  and 
observability  or  detectability  with  respect  to  extracting 


the  controlled  output  yc*  The  beam  state  and  the  target 
postlon  state  combine  to  form  the  controlled  output,  and 
detectability  can  be  shown  to  be  satisfied.  However,  as 
implied  in  Jamerson's  target  model,  there  is  assumed  to  be 
no  control  inputs  applied  (by  the  controller)  to  the  target. 
In  essence,  we  are  assuming  that  the  presence  of  the  beam  on 
the  target  has  no  impact  on  the  target  dynamics.  From 
linear  control  theory,  all  the  states  of  the  target  are  thus 
uncontrollable  from  the  point  of  entry  of  u  [3:444;  18:92). 
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in  addition,  the  presence  of  the  two  direct  Integrations  in 
the  dynamics  model  indicate  the  target  motion  of  the  open 
loop  system  is  astable.  This  is  the  same  as  having  two 
poles  at  the  origin  in  the  s-plane  or  two  poles  at  the  unity 
value  of  real  axis  of  the  z-plane.  Without  controllability 
or  stability,  there  is  no  guarantee  of  stabilizability  from 
the  point  of  entry  of  u,  nor  of  a  solution  to  the  steady- 
state  Riccati  equation  used  to  generate  constant  controller 
gains  (see  Section  3.2.2)  from  the  quadratic  cost  function. 
This  requires  a  small  adjustment  to  the  controller  model  for 
the  target.  By  moving  the  two  poles  at  the  origin  of  the 
s-plane  to  the  left  by  some  small  epsilon,  the  controller 
model  is  now  stable,  although  still  not  controllable.  Thus 
a  solution  to  the  steady-state  Riccati  equation  is 
guaranteed.  By  making  this  epsilon  small  and  at  least  an 
order  of  magnitude  away  from  the  time  constants  of  the 
target  and  the  beam,  the  dynamics  of  the  model  system  remain 
dominated  by  the  target  and  beam  time  constants.  Therefore 
the  gains  derived  from  the  solution  to  the  Riccati  equation 
should  still  be  applicable  to  the  filter  and  system  dynamics 
equations  involving  pure  integrations  C 8 ; 9 ] ,  though 
performance  evaluations  are  critical,  to  determine  the 
impact  of  this  model  alteration.  Thus,  the  adjusted 
controller  gain  model  is  given  by, 

Jt<ti+1)  *  4' <t1+1»ti)jL(t1)  +  ^(tjjuttj)  (3-58) 
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and  #'  is  determined  from  the  inverse  Laplace  transform  of 
Isl.-F*  1  and  the  adjusted  is  given  by: 
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(3-60) 


It  is  important  to  note  that  this  adjustment  is  only 
applicable  for  the  controller  model.  The  filter  and  truth 
models  remain  unchanged.  Vith  this  change  we  can  use  the 
techniques  of  Section  3.2.2  to  determine  the  appropriate 
gains . 
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3.3.3  Proportional  Plus  Integral  Regulator.  The 
proportional  gain  regulator  lacks  type  1  characteristics. 
Without  type  1  control  the  steady  state  (mean)  error  cannot 
be  driven  to  zero  when  constant  unmodelled  disturbances  are 
present.  To  achieve  type  1  control  we  need  to  include 
integral  action  in  the  control  loop.  Just  as  a  model  of  the 
beam  Is  incorporated  into  the  regulator  synthesis  model  so 
the  gains  are  appropriately  determined,  we  also  must 
Incorporate  a  model  of  the  Integration  process.  This  allows 
the  appropriate  feedback  gains  to  be  determined  for  each 
state  and  the  Integrator  by  LQ  synthesis  methodology,  rather 
than  a  trial  and  error  approach  to  determine  the  best  gain 


90 


settings.  Since  we  are  using  a  discrete  system,  it  is 
appropriate  to  model  the  integration  process  as  a  summation 
process  or  pseudointegral  (9:46;18).  Following  Johnson's 
development,  we  include  a  non-zero  setpoint  control  yd .  The 
nonzero  setpoint  control  will  allow  pointing  of  the  beam  at 
an  aimpoint  different  from  the  detector  array  centroid.  The 
pseudointegral  becomes  (18:133): 


qtv  *  q(t0)  +  z o  ly^tj)  -  ydi 


q(t,_1)  +  (xB(t1.1)  -  yd) 


( 3-61a ) 

( 3-61b) 


Figure  3.5a  illustrates  the  pseudointegral  for  the  PI 
regulator.  The  pseudointegral  is  included  in  the  controller 
model  by  augmenting  Equation  (3-61b)  to  controller  model  of 
Section  3.1.3,  Equation  (3-51).  The  augmented  determin¬ 
istic  controller  gain  model  becomes. 


% 

V141> 

* 

<v 

f  *  +  Vd 

( 3-62a ) 

i' 

01 

"x-it.n  [B.i  r< 

)' 

v,  ( 3-62b) 

•/ 

L«,tiu»  J 

c 

lj 

J  L°  J  1  Is 

L  ® 

The  optimal  control  law  is  given  by: 

u*(tj)  -  -£* (t1)jca(t1) 


(3-63) 


The  optimal  control  gain  is  calculated  as  is  Section  3.2.2. 
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3.3.4  Proportional  Plus  Integral  Tracker.  As  with  the 


PI  regulator,  we  also  Incorporate  a  model  of  the  integration 


process  in  the  PI  tracker.  This  allows  the  appropriate 


feedback  gains  to  be  determined  for  each  state  and  the 


integrator  by  LQ  synthesis  methodology,  rather  than  a  trial 


and  error.  As  before,  we  are  using  a  a  summation  process  or 


pseudointegral  as  a  model  of  the  integration  process.  Again 


we  include  a  non-zero  setpoint  control  yd .  However,  the 


nonzero  setpoint  control  will  now  allow  pointing  of  the  beam 


at  an  aimpoint  different  from  the  target  centroid.  The 


pseudointegral  then  becomes  t 9 : 46; 18 : 133 ] , 


q<V  *  q(t0)  +  lyB<V  -  x.j.pU.j)  -  ydl  ( 3-64a ) 


+  lxB(ti-l)  '  xTP(ti-l)  "  *d]  (3-64b) 


Figure  3.5b  illustrates  this  pseudointegral.  The  pseudo¬ 


integral  is  included  in  the  controller  model  by  augmenting 


Equation  (3-64b)  to  controller  model  of  Section  3.1.3 


Equation  (3-58).  The  augmented  deterministic  controller 


gain  model  now  becomes. 


VW  *  *aW  +  Bdau(ti>  +  Vd 


( 3-65a ) 
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yd  ( 3-65b) 


Equation  (3-63)  can  be  written  as  (with  time  arguments 
dropped) : 


(3-66) 


This  is  the  form  of  the  elemental  controller  used  in  this 
study.  Figure  3.6  depicts  this  controller  in  block  diagram 
form.  Elemental  controllers  of  this  form  will  be  used 
together  in  a  multiple  model  structure  described  next. 

3.3.5  Multiple  Model  Adaptive  Controller.  The  multiple 
model  adaptive  controller  (MMAC)  was  incorporated  by  Johnson 
to  combat  the  problem  sluggish  response  to  highly  dynamic 
targets.  The  derivation  of  the  MMAC  parallels  that 
presented  in  Section  2.4  for  multiple  model  adaptive 
estimation  Equations  (2-23)  through  (2-31).  Analogous  to 
that  development,  the  feedback  out  of  the  MMAC  is  the 
weighted  sum  of  the  outputs  of  each  elemental  controller  as 
given  by: 


subscript  "a"  denotes  "augmented".  For  Johnson's  MMAC  [9], 
each  p^  was  the  probability  weight  associated  with  each 


Kalman  filter.  However,  here  p^Ctj)  is  probability  of  both 
the  Heer  filter  estimate  and  the  Kalman  filter  estimate 
given  by: 

pk(ti}  “  pMeer(m)(tl)  pKalman( £) ( ti }  (3~69 


M 

pMeer (m)  =  1 

m=x 


(3-70) 


PKalman(i)  =  1  (3-71) 

where  here,  M  is  the  number  of  Meer  filters  and  here,  L  is 
the  number  of  Kalman  filters.  Thus,  K*LM.  Figure  3.7 
depicts  this  MMAC. 

3.3.6  Choosing  a  Controller.  In  previous  thesis 
efforts,  the  controllers  developed  were:  the  Proportional 
Gain  (PG)  controller  [32],  the  Proportional  plus  Integral 
(PI)  controller  (8;  23],  and  the  multiple  model  adaptive 
controller  (MMAC)  composed  of  a  bank  of  PI  controllers  [91. 
Each  of  these  controllers  were  constant  gain  devices.  To 
check  the  feasibility  of  using  the  state  estimates  provided 
by  the  Meer  filter  in  a  controller  design,  the  proportional 
gain  controller  was  developed  by  Zicker  (321.  This 
controller  featured  a  PG  target  tracker  designed  around  a  PG 
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beam  regulator.  VI th  this  design,  Zicker  shoved  that  the 
use  of  the  Meer  filter  in  a  controller  design  was  feasible. 
However,  the  PG  controller  inherently  lacks  the  type-one 
characteristics  necessary  to  reject  unmodelled  constant 
disturbances  injected  into  the  target  position  state  [3:196; 
18:132].  As  models  can  only  approximate  the  real  world,  our 
system  will  always  be  subject  to  these  unmodelled 
disturbances.  Thus,  Moose  [231  began  development  of  a 
controller  based  on  a  type-one  Proportional  plus  Integral 
(PI)  target  tracker  designed  around  a  beam  regulator. 
Jamerson  continued  Moose's  development,  but  exchanged 
Moose's  PI  tracker  with  another  PI  tracker  based  on  the  more 
realistic  target  model  that  Jamerson  incorporated  into  that 
study.  These  PI  trackers  were  able  to  reject  constant 
unmodelled  disturbances  (drive  the  mean  error  between  the 
beam  and  the  target  positons  to  zero),  as  desired.  However, 
they  were  slugglish  and  had  difficulty  tracking  highly 
maneuverable  targets  [8;  9).  To  improve  performance  against 
highly  maneuverable  targets,  Johnson  (9),  developed  a 
multiple  model  adaptive  controller  ( MMAC )  based  on  elemental 
controllers  that  used  a  PI  target  tracker.  Each  controller 
was  based  on  a  separate  assumption  about  the  approximate 
strength  for  the  assumed  target  dynamics  driving  noise. 
Although  variations  in  dynamic  driving  noise  strength  have 
no  effect  on  the  calculation  of  controller  gains,  the 
controller  outputs  differed  due  to  the  use  of  separate 
Kalman  filter  models  of  the  target  to  supply  the 
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corresponding  controllers  with  state  estimates  o£  the 
target.  Thus,  the  controller  gains  in  Figure  3.7  show  no 
dependency  on  the  individual  Q  value  assumed  by  the  target 
state  Kalman  filter.  The  differences  in  the  assumed  target 
dynamics  driving  noise  strength  does  affect  the  state 
estimates  generated  by  each  Kalman  filter.  Although  the 
MMAC  performed  well  [9],  Johnson  found  a  severe  robustness 
problem  in  the  elemental  PI  controllers  for  a  mismodelled 
beam  time  constant.  One  suggested  solution  to  the 
robustness  problem  is  the  use  of  multiple  model  adaptive 
techniques  on  the  beam  time  constant.  That  is  the  technique 
explored  in  this  study. 

Using  multiple  model  techniques  with  the  beam  time 
constant  results  in  a  structure  similar  to  Johnson's  MMAC . 
There  are  some  notable  differences  between  the  structures. 
Unlike  an  adaptation  on  the  target  driving  noise  strength, 
adapting  on  the  beam  time  constant  does  result  in  different 
controller  gains  for  each  controller.  Each  controller  then 
uses  the  beam  estimate  from  its  corresponding  Meer  filter, 
based  on  the  same  beam  time  constant,  to  generate  the 
appropriate  feedback.  Each  controller  output  is  then 
weighted  based  on  the  residual  performance  of  its  target 
filter  and  its  Meer  filter  to  form  the  actual  output.  Thus, 
the  controller  structure  for  this  study  is  a  bank  of  nine  PI 
controllers  based  on  a  unique  hypothesis  of  the  beam  time 
constant  and  target  driving  noise  strength.  As  in  past 
three  thesis  efforts,  each  elemental  PI  controller  is  based 


on  a  PI  tracker.  The  appropriate  constant  gains  are 
generated  through  the  use  of  LQ  synthesis  techniques  based 
on  the  assumption  of  full-state  feedback.  Then  through  the 
use  of  assumed  certainty  equivalence,  the  controller 
constant  gain  blocks  are  provided  the  target  and  beam 
estimates  from  the  respective  Kalman  and  Meer  filters.  Once 
we  have  chosen  a  structure  or  our  controller,  we  can  begin 
the  synthesis  of  the  controller  gains  via  LQG  and  Implicit 
model  following  techniques.  Section  3.2  in  fact  discussed 
this  synthesis  in  detail. 

3.4  Summary 

This  chapter  presented  the  current  target  model  being 
used  in  this  research  as  a  third  order  Gauss-Markov  position 
process.  This  satisfied  the  linear  model  and  Gaussian  noise 
requirements  for  a  tractable  control  design.  Next,  the 
quadratic  cost  function  for  synthesizing  each  elemental 
controller  gain  was  described.  The  quadratic  cost  function 
results  in  a  tractable  approach  to  controller  design. 

Linear  models  and  a  quadratic  cost  criterion,  allows  the  use 
of  LQ  synthesis  techniques  to  determine  controller  gains. 
Assumed  certainty  equivalence  is  then  used  to  incorporate 
these  gains  into  the  NPB  controller.  Weighting  matrices  in 
the  cost  function  allow  the  user  to  specify  how  the  control 
energy  will  be  applied  to  regulate  the  states.  Implicit 
model  following  techniques  are  one  method  of  determining  the 
appropriate  weighting  matrices  for  acheiving  desired 
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performance  and  robustness.  The  controller  structure  chosen 
for  this  study  was  a  multiple  model  PI  controller.  The  PI 
controller  drives  steady  state  mean  errors  to  zero  even  in 
the  presence  of  constant  unmodelled  disturbances.  The 
multiple  model  approach  attacks  performance  and  robustness 
problems  by  processing  a  number  of  simple  controllers  in 
parallel,  rather  than  increasing  the  complexity  of  a  single 
controller.  The  next  chapter  discusses  the  analysis 
performed  on  the  controller  design  of  this  chapter. 


IV 


.  SIMULATION  DESIGN  AND  EVALUATION  PROCEDURE 

A  design  is  not  complete  until  its  performance  is 
validated  against  some  pre-specif ied  standard,  baseline, 
goal,  or  specification.  The  first  step  in  this  validation 
process  is  an  evaluation  of  performance  using  computer 
simulation.  This  chapter  presents  the  method  used  to 
simulate  the  NPB  tracking  problem  and  evaluate  estimator  and 
controller  performances  in  this  study. 

Two  performance  evaluations  are  made.  The  first 
evaluation  is  of  the  MMAE  Meer  filter's  ability  to  estimate 
the  position  of  the  NPB  centroid.  The  MMAE  Meer  filter  is 
composed  of  a  bank  of  Meer  filters.  Each  Meer  filter  is 
based  on  a  separate  hypothesis  about  the  value  of  the  beam 
time  constant.  The  second  evaluation  is  of  the  ability  of 
an  MMAC  based  on  an  implicit  model  following  to  point  the 
NPB  at  a  given  target.  The  MMAC  is  composed  of  a  bank  of 
nine  PI  controllers,  each  operating  with  a  unique 
combination  of  one  of  three  Kalman  filters  and  one  of  three 
Meer  filters.  Each  Kalman  filter  is  based  on  a  separate 
hypothesis  about  the  strength  of  the  target  driving  noise. 

This  chapter  begins  by  presenting  the  analysis 
procedure.  Following  this,  it  describes  the  software  tools 
used  to  produce  and  assimilate  the  data  for  performance 
evaluations.  Finally,  filter,  controller,  and  simulation 
parameters  used  in  this  study  are  presented. 
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4.1 


4.1  Performance  Analyses 

Eight  different  performance  analyses  are  required  to 
evaluate  the  MMAE  Meer  filter  and  to  evaluate  the  use  of 
implicit  model  following  weighting  matrices  in  the  synthesis 
of  controllers.  The  first  five  analyses  involve  the 
evaluation  of  the  MMAE  Meer  filter.  The  last  three  analyses 
involve  the  MMAC  using  cost  weighting  matrices  developed  by 
Implicit  model  following  techniques. 

The  MMAE  Meer  filter  is  an  adaptive  filter  composed  of  a 
bank  of  elemental  Meer  filters,  each  based  on  an  assumed 
beam  time  constant  7b.  The  Meer  filter  itself  is  an 
adaptive  filter  composed  of  elemental  Synder -Fishman 
filters,  each  based  on  a  different  hypothesis  about  the 
sequence  of  noise  and  signal  events.  The  effects  of  such  a 
second  level  of  adaptation  are  unknown.  Therefore,  an 
evaluation  of  the  estimation  and  control  performance  with 
the  MMAE  Meer  filter  is  necessary  to  determine  whether  this 
adaptation  can  adequately  estimate  the  beam  location  for 
off-nominal  beam  time  constants. 

The  MMAC  is  composed  of  a  bank  of  elemental  PI 
controllers.  Each  elemental  controller  uses  one  of  three 
Kalman  filters  to  estimate  the  target  state  (differing  from 
each  other  on  basis  of  dynamics  driving  noise  strength)  and 
one  of  three  Meer  filters  (differing  on  assumed  beam  time 
constant)  to  estimate  the  beam  state.  Bach  PI  controller 
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will  be  based  on  a  unique  combination  o£  the  target  dynamics 
driving  noise  strength  [9]  and  beam  time  constant.  Thus  the 
bank  will  be  composed  of  9  elemental  PI  controllers.  The 
analysis  of  the  controller  performance  will  allow  selection 
of  the  weighting  matrices  (used  for  controller  synthesis) 
that  provide  the  best  tradeoff  of  performance  at  design 
conditions  versus  robustness. 

4.1.1  Evaluating  the  MMAE  Meer  Filter.  Five  analyses 
will  be  conducted  to  evaluate  the  performance  of  the  MMAE 
Meer  filter.  The  first  analysis  set  will  determine  whether 
the  use  of  an  MMAE  Meer  filter  can  consistently  provide 
precise  estimates  of  the  beam  location  when  the  beam  time 
constant  is  uncertain.  That  a  MMAE  Meer  filter  provides 
good  state  estimates  does  not  imply  that  it  also  provides 
good  parameter  estimates.  The  second  analysis  will 
determine  whether  the  MMAE  Meer  filter  can  provide  an 
adaptive  estimate  of  the  beam  time  constant  while  providing 
good  state  estimates.  As  such,  this  analysis  is  run 
concurrently  with  the  first.  The  third  analysis  will 
justify  the  use  of  an  elemental  Meer  filter  depth  of  one 
(D  =  1 ) ,  The  last  two  analyses  will  define  the  limits  on 
the  beam  time  constant  for  which  the  MMAE  Meer  filter  will 
be  able  to  operate. 

Johnson  experienced  significant  degradation  in  target 
tracking  performance  with  only  2%  variation  in  filter  beam 
time  constant  from  the  true  value  [9:1181.  Since  he  found 


little  change  In  the  controller  gains  for  small  changes  In 
the  beam  time  constant  (9:124],  we  suspect  that  the  Meer 
filter  does  not  provide  an  adequate  estimation  of  the 
position  of  the  beam  centroid  position  when  the  beam  time 
constant  is  mismodeled.  Thus  the  focal  point  for  our  first 
analysis  is  the  MMAE  Meer's  ability  to  estimate  NPB  centroid 
location. 

For  the  first  analysis,  the  MMAE  Meer  filter  consists  of 
3  Meer  filters,  based  on  the  nominal  beam  time  constant, 
nominal  minus  4%,  and  nominal  plus  4\,  with  a  Meer  filter 
memory  depth  of  3.  The  weighted  beam  position  estimate  of 
the  filter  is  compared  to  the  true  beam  position  to  evaluate 
true  performance  capability.  This  first  analysis 
demonstrates  whether  a  simple  MMAE  Meer  filter  can 
consistently  provide  accurate  beam  position  estimates  when 
the  true  beam  filter  time  constant  differs  from  the  assumed 
time  constants  of  the  individual  models.  The  true  value  of 
the  beam  time  constant  will  be  set  to  the  nominal,  at 
nominal  plus  2%,  and  at  nominal  minus  2%.  The  first  setting 
provides  an  optimal  setting  of  the  true  beam  time  constant 
with  respect  to  the  filter-assumed  values.  The  last  two 
setting  are  to  provide  a  worst  case  performance  by 
mismatching  the  true  value  with  the  filter-assumed  values  as 
much  as  possible:  no  single  elemental  filter  is  specifically 
tuned  for  these  parameter  values.  If  performance  is  poor, 
then  additional  Meer  filters  may  be  required  to  provide  a 
finer  discretization  of  the  beam  time  constant.  This 


analysis  will  indicate  whether  the  MMAE  Meer  filter  is  a 
viable  method  of  solving  the  robustness  problem  found  by 
Johnson  with  regard  to  the  beam  time  constant  [9]. 

Subsequent  analyses  will  use  an  MMAE  with  the  number  of  Meer 
filters  determined  as  appropriate  in  this  first  analysis. 

The  second  analysis  determines  whether  the  filter  can 
provide  an  adaptive  estimate  of  the  beam  time  constant  while 
providing  good  state  estimates.  A  good  adaptive  beam  time 
constant  estimate  could  be  used  to  propagate  the  beam 
estimate.  A  good  estimate  of  the  beam  could  also  be  used  in 
adaptive  controller  designs  that  calculate  the  feedback  gain 
in  real  time  based  on  an  estimate  of  the  beam  time  constant. 

The  third  analysis  evaluates  the  performance  of  the  MMAE 
Meer  filter  when  the  filter  memory  depth  is  decreased  from 
D  *  3  to  D  *  1.  This  will  reduce  the  number  of  propagating 
elemental  Snyder -Fishman  filters  in  each  elemental  Meer 
filter  from  eight  to  two.  Additionally,  each  elemental  Meer 
filter  can  be  expressed  equivalently  with  only  a  single 
elemental  Synder -Fishman  filter,  but  with  its  gain  expressed 
as  a  function  of  residual  size,  as  seen  in  Section  2. 5. 2. 2 
(6;  9;  32].  Although  this  may  slightly  increase  the  average 
RMS  error,  it  should  significantly  decrease  the 
computational  loading.  It  makes  on-line  feasibility  much 
more  reasonable,  and  so  a  tradeoff  analysis  is  warranted. 

The  fourth  MMAE  Meer  filter  analysis  will  determine  the 
effects  of  a  linearly  changing  and  sinusoidally  varying  true 
beam  time  constant  on  the  performance  of  the  MMAE  Meer 


filter.  This  will  give  insight  into  the  MMAE  Meer  filter's 
ability  to  track  a  real  world  particle  beam  with  physically 
motivated  parameter  changes. 


The  final  analysis  will  determine  the  limits  on  the 
variation  of  beam  time  constant  using  the  MMAE  Meer  filter. 
Limits  using  additional  Meer  filters  may  be  investigated. 
These  limits  would  Indicate  the  discretization  necessary  to 
cover  the  range  of  beam  time  constants  over  which  the  MMAE 
Meer  filter  may  operate. 

4.1.2  Evaluating  the  MMAC  with  Implicit  Model  Following 
Weighting  Matrices.  In  the  first  controller  analysis, 
candidate  weighting  matrices  are  determined  using  Implicit 
model  following.  The  performance  with  these  weighting 
matrices  is  measured  to  ensure  reasonable  performance  at 
design  conditions.  The  next  two  analyses  indicate 
performance  for  variations  in  beam  and  target  parameters 
with  the  weighting  matrices  that  have  reasonable  performance 
at  design  conditions  as  determined  in  the  first  analysis.  A 
robustness  analysis  is  performed  to  indicate  performance 
when  controller-assumed  parameters  do  not  match  actual  beam 
and  target  parameter.  This  analysis  will  allow  the 
selection  of  the  weighting  matrices  that  provide  the  best 
tradeoff  in  performance  at  design  conditions  versus 
robustness.  A  sensitivity  analysis  is  performed  to  indicate 
performance  capabilities  when  the  controller  is  operated 
over  a  variety  of  design  conditions.  In  a  sensitivity 


analysis,  the  controller-assumed  parameter  matches  the 
corresponding  beam  or  target  parameter.  This  will  indicate 
the  range  of  design  conditions  over  which  the  weighting 
matrices  are  applicable. 

The  determination  of  the  best  cost  weighting  matrices 
based  on  implicit  model  following  techniques  is  not  an  exact 
science.  There  may  exist  a  number  of  sets  of  weighting 
matrices  which  provide  an  adequate  level  of  performance  with 
a  reasonable  degree  of  robustness.  Therefore,  a  performance 
analysis  is  conducted  to  determine  the  weighting  matrix  set 
which  provides  good  performance  at  design  conditions  and 
also  enhances  the  MMAC's  robustness.  In  other  words,  the 
goal  is  to  select  the  set(s)  of  weighting  matrices  that  best 
minimize  the  average  RMS  error  (without  using  excessive 
amounts  of  control)  when  the  truth  model  matches  the  filter 
model,  while  simultaneously  providing  at  least  stability  and 
some  degree  of  desirable  performance  when  controller-assumed 
parameters  do  not  match  corresponding  beam  and  target 
parameters . 

A  robustness  analysis  is  performed  to  evaluate  the  MMAC 
controller's  performance  when  a  controller-assumed  model 
parameter  differs  from  that  of  the  truth  model.  The 
robustness  analysis  is  conducted  under  the  same  guidelines 
as  the  sensitivity  analysis,  but  differs  in  that  a  parameter 
is  changed  in  the  truth  model  without  informing  the 
controller.  The  purpose  of  this  study  is  to  determine 
whether  the  use  of  implicit  model  following  techniques  in 
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controller  design,  can  result  In  the  design  of  a  robust 
controller  while  maintaining  adequate  performance  at  design 
conditions . 

The  baseline  against  which  these  analyses  are  compared 
is  defined  as  a  PI  controller  which  receives  the  beam  and 
target  states  from  the  truth  model  and  has  access  to  the 
truth  model  parameters.  Such  a  baseline  provides  the  best 
performance  possible  for  the  chosen  set  of  controller  gains. 
A  second  benchmark  is  defined  as  a  PI  controller  which 
receives  state  estimates  from  a  single  Kalman  filter  and  a 
single  Meer  filter  which  have  access  to  the  correct 
parameters.  The  second  baseline  provides  a  realistic 
assessment  of  the  best  performance  possible  from  controllers 
that  operate  with  state  estimates  rather  than  the  actual 
states . 

The  purpose  of  the  sensitivity  analysis  is  to  evaluate 
the  performance  of  the  MMAC  controller  based  on  the 
weighting  matrices  determined  in  the  previous  analysis  as 
its  design  parameters  (both  truth  and  filter)  are  changed. 
Sensitivity  to  parameter  variations  may  indicate  the  need 
for  further  adaptive  estimation.  Since  the  contoller  gains 
are  dependent  on  the  design  parameters,  as  well  as  choice  of 
weighting  matrices,  robust  performance  at  a  given  set  of 
design  parameters  may  not  be  indicative  of  robustness  for 
other  sets  of  design  parameters. 

The  sensitivity  analysis  is  based  on  the  variation  of  a 
single  parameter  setting  in  both  the  controller-assumed 
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model  and  truth  model  (i.e.  "real  world"),  while  other 
parameters  remain  at  the  nominal  values  presented  in  the 
previous  sections.  For  this  study,  data  will  be  collected 
as  parameters  are  varied  one  order  of  magnitude  above  and 
below  the  nominal  value.  The  analysis  evaluates  all  the 
beam  and  target  parameters  except  the  filter  depth,  which  is 
set  to  one  (D  *  1).  These  results  are  compared  to  the 
results  obtained  by  Johnson  [9:951  using  the  cost  weighting 
matrices  that  are  based  on  the  results  in  Chapter  5  of 
Johnson's  thesis.  These  weights  are:  U=l,  Xii=100,  X55=10 
and  all  other  weights  and  cross  weight  equal  to  zero,  where 
U  is  the  weight  on  the  control  energy,  Xn  is  the  weight  on 
tracking  error  and  X55  is  the  weight  on  the  pseudo¬ 
integration  state  (See  Section  3.3.4). 

4.1.3  The  Method  -  Monte  Carlo  Simulation.  It  will  be 
necessary  to  evaluate  the  MMAE  Kalman  filter  and  MMAE  Meer 
filter  state  estimation  error  statistics,  to  ensure  these 
filters  are  not  diverging  and  are  in  fact  providing  viable 
estimates  of  important  state  variables.  It  will  also  be 
required  to  evaluate  MMAC  tracking  errors.  The  performance 
of  an  estimator/controller  can  often  be  analyzed  by  either 
covariance  analysis  or  Monte  Carlo  simulation  [16:3291.  The 
least  time  consuming  method  is  covariance  analysis,  which 
requires  only  one  software  run  to  generate  the  time  history 
of  the  estimation  error  covariance.  Unfortunately, 
covariance  analysis  should  not  be  applied  to  this  system,  as 


proper  covariance  analysis  is  limited  to  linear  stochastic 
estimators  and/or  controllers  that  use  prespecified 
measurement  update  times.  The  adaptive  mechanisms  in  both 
the  controller  and  the  MMAE  Meer  filter,  as  veil  as  the  Meer 
filter  itself,  violate  the  linearity  restriction,  and  the 
elemental  Synder -Fishman  filters  have  a  varying  sample  rate 
that  is  not  prespecif iable .  Therefore,  the  more  time 
consuming  Monte  Carlo  simulation  is  used  to  evaluate  the 
particle  beam  estimator/controller.  Next,  we  describe  the 
software  used  to  facilitate  the  generation  and  collection  of 
data  for  this  study. 


The  Tools  -  SOFE  and  SOFEPL 


The  performance  data  are  generated  from  a  specially 
modified  version  of  a  software  package  called  SOFE 
(Simulation  for  Optimal  Filter  Evaluation)  [24].  Subsequent 
performance  analysis  is  performed  using  SOFEPL  (SOFE 
Plotter)  (61.  SOFE  is  a  general  purpose  Monte  Carlo 
simulation  program  designed  for  evaluating  systems  that  use 
Kalman  or  extended  Kalman  filters.  This  section  begins  by 
describing  how  SOFE  was  modified  and  used  to  simulate  the 
NPB,  target,  Meer  filter,  Kalman  target  filter  and  the 
controller.  Following  this,  the  use  of  SOFEPL  to  assimilate 
the  data  into  error  statistics  is  presented. 


SOFE  allows  the  user 


to  simulate  a  truth  model  and  a  Kalman  or  an  extended  Kalman 


filter  model  that  can  be  described  by  a  set  of  stochastic 
differential  equations.  SOFE  also  allows  the  user  to 
specify  the  measurement  format  of  the  filter  update,  and  if 
desired,  to  include  impulsive  feedback  control.  However, 
the  incorporation  of  the  Meer  filter  required  a  modification 
to  the  basic  SOFE  code.  In  Kalman  or  extended  Kalman  filter 
applications,  the  update  times  occur  at  prespecif iable 
constant  intervals,  whereas  the  Meer  filter  updates  occur  at 
random  intervals.  Therefore,  it  was  necessary  for  Meer  to 
modify  the  SOFE  software  to  allow  propagation  to  the  Meer 
filter  update  time.  Therefore,  this  section  will  first 
describe  the  modifications  necessary  to  include  random 
Interval  updates.  Then,  the  equations  used  to  simulate  the 
NPB  and  target  models  are  presented.  Following  this,  the 
routines  necessary  to  simulate  the  Meer  filter  model  are 
described.  Then  the  Kalman  filter  model  for  the  target  is 
presented.  Finally,  the  method  used  to  provide  continuous 
feedback  (piecewise  constant,  as  generated  by  a  computer  and 
passed  through  a  zero  order  hold)  with  the  MMAC  controller 
is  described.  The  SOFE  source  code  that  Meer  altered  and 
that  is  used  throughout  the  research  is  the  FORTRAN  4,  May 
1982  version  of  SOFE  (21:791. 

4. 2. 1.1  Random  Interval  Updating.  Since  the 
Kalman  filter  is  based  on  a  fixed  sample  period,  SOFE 
updates  the  Kalman  filter  at  regular  time  intervals.  Thus, 
Meer  (211  had  to  modify  the  SOFE  source  code  to  accommodate 


» 


112 


the  random  arrival  time  of  Poisson  distributed  signal  and 
noise  events  at  the  detector  (See  Section  2.2).  Meer 
altered  the  time  that  a  user-defined  error  statistics 
generation  routine  is  called,  from  a  constant  interval  to 
the  computed  time  of  the  next  signal  or  noise  event, 
whichever  is  soonest  to  occur.  For  simplicity,  the  times  to 
the  next  signal  and  noise  events  are  calculated  and  compared 
in  the  simulation  from  user-defined  signal  and  noise  rate 
parameters.  The  equation  used  to  calculate  the  varying 
sample  period  is  derived  from  the  Poisson  process  density 
function. 


p(y)  =  -([Atly/yl]  expt-Atl  (4-1) 

where  p(y)  is  the  probability  that  y  events  occurred  within 
time  period  t,  and  A  is  the  mean  signal  arrival  rate  or  mean 
noise  arrival  rate.  Because  we  are  interested  in  generating 
a  random  sample  period  with  the  next  event  occurring  at  the 
completion  of  the  sample  period,  y  is  set  to  zero  (i.e.,  no 
events  prior  to  the  end  of  the  sample  period),  and  the 
random  sample  period,  t,  is  calculated  from  the  inverse 
mapping  function  of  Equation  (4-1)  (2;  9:83-84;  29:621. 

t  =  -lAl"1  ln(p(0) 1  (4-2) 


where  p(0)  Is  defined  as  the  probability  that  t  amount  of 
time  has  passed  before  an  event  arrives,  and  01*1.  This 
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probability  has  a  range  o£  0<p(0)<l.  Because  we  are 
interested  in  generating  samples  of  tj,  p(0)  is  selected 
randomly  from  the  range  of  p(0),  and  this  process  is  used  to 
simulate  the  arrival  time  of  the  next  signal-  or 
noise-induced  event.  Both  the  signal  and  noise  are  Poisson 
processes ,  but  each  has  its  own  mean  arrival  rate  and  they 
are  Independent  processes.  Thus,  modified  SOFE  simulates  a 
varying  sample  period  (the  time  to  the  next  signal-  or 
noise-induced  event,  whichever  comes  first)  as  a  Poisson 
time  process. 

The  mean  signal  arrival  rate  is  defined  as 


As  -  r  ( 2*R)1/2  =  fi/(tF-tQ) 


(4-3) 


where  r  is  the  maximum  amplitude  of  the  signal  rate  density 
function,  R  is  the  beam  dispersion,  and  B  is  the  number  of 
signal-induced  events  expected  over  the  duration  of  the 
simulation  run,  (tp  -  tg).  The  relationship  between  the 
mean  signal  arrival  rate  and  the  expected  signal  rate 
parameter,  A_(t,r,x(t) ),  can  be  shown  by  integrating 
Equation  (2-7)  (shown  for  the  one  dimensional  case). 


Ag(t1f«/X(t) )d« 


r(t^)exp(-aTR_1(t1)«/2Ida 
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where  a=tr-H(t)x(t ) ] .  For  this  study,  we  are  assuming  that 
r  and  R  are  constant.  Although  there  Is  no  known  closed 
form  expression  for  the  integral  of  the  Gaussian  density 
function,  we  know  that  the  entire  area  under  a  probability 
density  function  must  have  a  magnitude  of  one. 


J  t  ( 2n  r  r1/2 


exp(-«T  R-1  «/2)dot 
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(4-5) 


where  m  is  the  dimensionality  of  the  vector  a  and  equals  one 
for  this  study.  Although  our  detector  does  not  extend  to 
iniflnity,  we  can  assume  that  the  area  lost  is  neglible  if 
the  detector  width  is  much  larger  than  the  beamwidth  and  the 
beam  is  kept  away  from  the  edges  of  the  detector.  Thus,  the 
maximum  amplitude  of  the  rate  function,  r,  is: 


r  =  ia/(tF-t0)]  U2*)  R  ]~1/2 


(4-6) 


The  mean  noise  arrival  rate  is  defined  as 


(t,r )dr 


(4-7) 


where  L  is  the  physical  length  of  the  detector.  For  this 
study  An(t,r)  is  assumed  to  be  a  constant  denoted  as  An,  so 
that  X p ®  An  L. 


4. 2. 1.2  Truth  Models.  Adapting  SOFE  to  this 


application  did  not  affect  the  propagation  equations  for  the 
SOFE  truth  models  used  to  simulate  the  NPB  and  target.  In 


this  study  the  truth  model  propagation  equations  are  linear 
and  given  by  Equations  (2-3)  and  (3-4).  SOFE  uses  a  fifth 
order  Kutta-Merson  integrator  to  solve  the  deterministic 
portion  of  the  propagation  equations.  Then  a  stochastic 
term  is  added  to  user-specified  states  at  periodic  intervals 
(after  the  completion  of  each  integration  to  that  time) 
through  separate  calls  to  a  Gauss  random  number  generator 
for  each  state.  Expressed  in  equation  form: 

2tg  =  JJq  +  GAUSS  ( 0.,  Q^)  (4-8) 

where  xg  is  the  stochastic  solution,  x  D  is  the  deterministic 
solution,  and  the  expression  GAUSS (0.,  Q.^)  is  a  vector  of 
randomly  generated  noises  of  mean  zero  and  a  covariance  of 
ad-  The  discrete-time  dynamics  driving  noise  covariance, 
(the  second  moment  of  the  equivalent  discrete  time  noise 
representing  the  effects  of  G(t)j£(t)  on  a  state  over  a  user 
specified  propagation  interval),  is  calculated  by  the  user 
from  Equation  (2-19)  in  Section  2.3  for  the  beam  and  by 
Equation  (3-21)  in  Section  3.2.1  for  the  target. 

4. 2. 1.3  Target  Kalman  Filter.  SOFE  propagates  the 
the  Kalman  filter  states  in  to  same  manner  as  the  truth 


states  (see  previous  section).  The  equation  used  to 
propagate  the  filter  covariance  £.  in  SOFE  is  assumed  to  be 
of  the  form  (24:83): 


As  with  the  state  propagation  equations,  SOFE  uses  a 
fifth-order  Kutta-Merson  integrator  to  solve  Equation  (4-9) 
SOFE  uses  a  recursive  method  to  incorporate  the  Kalman 
filter  measure  updates.  In  this  study  the  measurement  model 
is  also  linear  and  given  by  equation  (3-8).  SOFE  injects 
the  measurement  noise  into  each  measurement  by  the  following 
relationship: 

z(t1)  =  GAUSS(H2Lrp(ti  ),  R)  (4-10) 

where  R  is  the  variance  of  the  target  sensor  measurement 
noise.  For  this  problem  the  measurement  noise  is  assumed  to 
be  mean  zero.  The  standard  Kalman  filter  updates  at  regular 
discrete  time  intervals  determined  by  the  user,  whereas  the 
Meer  filter  updates  only  when  it  receives  a  Poisson  point 
process  event. 

4. 2. 1.4  Meer  Filter.  Unlike  the  method  SOFE  uses 
to  propagate  the  (possibly  nonlinear)  Kalman  filter  states, 
the  Synder -Fishman  filters  which  make  up  the  Meer  filters 
are  propagated  using  the  state  transition  matrix.  This  is 
possible  only  because  of  the  simple  beam  model  we  have 
chosen  (see  Chapter  2).  The  propagation  is  given  by  (again 
in  one  dimension): 

x(tj+i)  -  ♦  x(tj)  +  Bd  u  (4-11) 

where  ♦  ■  exp(-dt/TB)  and  Bd  *  rB  (1-*)B  as  given  by 
Equation  (3-19),  and  B  is  constant  for  this  study.  The  time 


Interval  t  is  the  time  since  the  last  Poisson  process  event 
occurred.  The  Integrator  for  the  truth  model  and  the  Kalman 
filter  model  are  also  provided  the  time  of  the  next  Meer 
filter  update.  However,  the  propagation  of  the  truth 
states,  the  Kalman  filter  states  and  the  filter-computed 
error  covariances  is  accomplished  by  integrating  the 
differential  equations  from  the  current  time  forward  to  some 
specified  time.  This  time  can  be  the  integration  step  size, 
the  time  for  the  next  Snyder-Flshman  filter  update,  or  the 
time  for  the  next  Kalman  filter  update.  Thus,  the  Meer 
filters  propagate  to  all  the  predetermined  points  for 
control  application  and  also  to  all  the  measurement  event 
times.  The  Meer  filter  estimates  need  to  be  taken  at 
regular  Intervals  for  ensemble  processing.  Therefore,  they 
are  recorded  at  the  Kalman  filter  update  times. 

As  mentioned  before,  the  MMAE  Meer  filter's  measurement 
update  occurs  whenever  an  event  is  observed.  If  the  event 
was  signal-induced,  the  true  location  of  event  and  the 
measurement  is  given  by: 

z(t)  -  GAUSS ( H&pg ( t ) ,  R)  (4-12) 

where  HxTB  is  the  true  location  of  the  beam  centerline  on 
the  detector  and  v/r~  is  the  beam  halfwidth.  Ve  have  assumed 
that  the  detector  provides  perfect  measurements  in  this 
study.  If  the  event  is  noise-induced.  Equation  (4-12)  does 
not  apply  and  the  location  of  the  noise-induced  event  is 


simulated  by  using  a  uniformly  distributed  mapping  function, 
uniformly  distributed  over  the  length  L  of  the  detector. 


The  updating  to  the  Meer  Filter  and  the  MMAE  Meer  filter 
follows  the  method  given  in  Section  2.3.  This  is  followed 
by  the  pruning  of  the  Meer  filter  as  described  in  Section 
2.5.2.  The  last  major  portion  of  user  modification  is  that 
related  to  the  setup  and  operation  of  the  controller. 

4. 2. 1.5  Controller  Simulation.  Zicker  [321  added 
the  first  controller  design  as  discussed  in  Section  3.3.6. 
Moose  and  Jamerson  (8;  23]  designed  the  follow-on 
proportional  plus  integral  controllers  and  Johnson 
incorporated  the  first  MMAC  design  employing  proportional 
plus  integral  controllers.  All  of  these  designs  used  the 
basic  structure  developed  by  Zicker  [32:391  for  calling  the 
controller  routine.  As  such,  it  operates  independently  of 
the  MMAE  Meer  filter  propagation  and  update  cycle  and  the 
Kalman  filter  propagation  and  update  cycle  to  provide 
control  feedback  at  user-specified  intervals  using  the 
estimates  of  each  filter.  Normally  the  controller  update 
cycle  is  tied  to  the  Kalman  filter  cycle.  This  same  calling 
structure  is  used  in  this  study  so  that  results  can  be 
compared  to  Johnson's  work. 

4. 2. 1.6  Additional  Alterations.  The  only  other 
significant  modification  to  the  basic  SOFE  code  is  Zicker's 
alteration  of  a  portion  of  the  Gauss  random  number  generator 
to  increase  efficiency  (32:42).  The  conceptual  operation  of 
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the  modified  SOFE  program  is  shown  at  a  macro  level  In 
Figure  4-1.  SOFE  Is  used  to  generate  the  data  necessary  to 
determine  performance;  however  it  would  be  impractical,  if 
not  impossible,  to  determine  the  actual  performance  from  the 
raw  data.  Therefore,  another  software  package  is  used  to 
post-process  the  data. 

4.2.2  Assimilating  the  Data  -  SOFEPL.  SOFEPL  is  a 
postprocessing  routine  that  can  be  programmed  to  perform 
various  statistical  functions,  such  as  calculate  the  RMS 
error  time  history,  and  then  plot  the  results.  The  data 
provided  by  SOFE  is  statistically  reduced  and  the  results 
are  plotted  by  SOFEPL.  First  this  section  will  discuss  the 
statistics  used  to  evaluate  the  performance  of  the  MMAE  Meer 
filters.  Although  the  SOFEPL  program  provides  an  array  of 
statistical  options,  they  are  limited  to  the  error 
statistics  found  between  the  truth  and  filter  states  (for 
filter  performance  evaluation)  and  are  inadequate  for 
evaluating  a  controller  design.  An  alternative  approach 
provided  by  SOFEPL  is  for  the  user  to  specify  the  desired 
error  to  be  evaluated  statistically.  The  final  portion  of 


this  section  discusses  the  procedure  used  to  generate  the 
controller  RMS  error  time  history,  where  the  statistics  are 
based  on  the  error  between  the  true  beam  and  true  target 
positions . 


4. 2. 2.1  MMAE  Meer  Filter  Statistics.  The  MMAE 
Meer  filter's  state  estimation  error  statistics  are  used  to 
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filter  position  estimate  defined  in  equation  (2-15)  and  n  is 
the  Monte  Carlo  run  number.  Thus,  the  mean  error  time 
history  is: 


n»e  ( t  j  )  =  (1/N)  E  efi(ti,n) 
n=l 


(4-14) 


where  N  is  the  total  number  of  runs.  The  variance  and  the 
standard  deviation  of  the  error  are  [4:74-77;  6]: 


ve(ti)  =  ( 1/ ( N-l ) ] [ E  eg( tj,n) J  -  IN/(N-1)1  mgttj)  (4-15) 

n=l 


ae(ti)  =  (1  +  1/1 4 (N-l)  1 1  lve(t1)J 


(4-16) 


Equation  (4-15)  calculates  the  unbiased  estimate  of  the 
variance  for  all  the  runs  and  Equation  (4-16)  is  an 
approximation  that  produces  an  unbiased  (to  first  order) 
estimate  of  the  standard  deviation  for  all  the  runs.  A  full 
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derivation  of  the  exact  equation  and  the  limitations  that 
apply  to  the  approximation  can  be  found  in  Demlng  [5;  9:751. 

The  root-mean-square  (RMS)  tracking  error,  RMS^t^),  and 
the  time-averaged  RMS  tracking  error,  RMSQ,  are  used  to 
further  compress  the  statistics  used  to  evaluate  the 
performance  of  the  estimater/controller  during  the  different 
analyses.  The  time  averaged  RMS  error  (RMSe)  is  useful  as 
it  provides  a  single  number  indicating  both  the  precision 
and  consistency  of  estimator/controller  performance.  This 
single  number  simplifies  the  task  of  comparing  sets  of 
statistics.  However,  the  mean  error  plus  and  minus  one 
standard  deviation  time  histories  are  really  the  most  useful 
Information  for  evaluation  of  the  estimator/controller. 

The  RMS  error  time  history  can  be  generated  by  using  the 
equation: 

RMSe(ti)  =  (m^tj)  +cr^(ti)l1/2  (4-17) 

The  time  averaged  RMS  error  is  calculated  by  averaging  each 

RMS  error,  RMS  (t.),  from  an  initial  time  t  to  a  final  time 
el  0 

tF  by  the  equation: 

f 

RMS.  =  C 1/C f +  1 ) I  E  RMS  ( t . )  (4-18) 

e  i-0  e  1 

Similar  to  the  MMAE  Meer  filter  error  statistics  arc  the 
controller  error  statistics  discussed  next. 


tracker  problem  Is  defined  as  the  difference  between  the 
true  target  state  and  the  true  beam  position  state: 

eT(tj,n)  =  xtXp(ti'n)  ~  xtB(ti'n)  (4-19) 

where  x^Xp  is  true  target  position  generated  by  the  target 

truth  model  defined  in  Section  3.1  and  is  true  beam 

position  generated  by  the  NPB  truth  model  defined  in  Section 

2.1.  The  errors  from  each  run  are  averaged  and  a  mean  error 

time  history,  m^ft^),  a  variance  time  history,  ve(t^),  a 

standard  deviation  of  error  time  history,  or.(t.),  the  RMS 

*  1 

error  time  history,  and  the  time  averaged  RMS  error  are 
calculated  using  Equations  (4-13)  through  (4-19)  as  with  the 
beam  error.  Equation  (4-19)  is  also  used  to  evaluate  the 
filter  computed  tracking  error  by  replacing  the  actual 
target  and  beam  states  with  the  estimated  states  from  the 
MMAE  Meer  and  Kalman  filters.  The  actual  rms  error  that 
occurs  is  dependent  on  the  parameter  values  used  in  the 
truth  models,  filters  and  controller.  The  next  section  will 
present  the  parameters  used  in  this  study. 


Filter.  Controller,  and  Simulation  Parameters 


The  particle  beam  estimator/controller  problem  is 
modeled  as  a  stochastic  process  and  the  performance  can  be 
evaluated  by  observing  the  statistical  behavior  of  the 
estimation  and  control  tracking  error  processes.  We  are 


interested  in  both  the  sensitivity  and  the  robustness  of  our 
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estimator  or  controller  design.  In  evaluating  the 
controller's  sensitivity  and  robustness,  we  would  like  to 
analyze  tracking  error  statistics  for  all  physically 
reasonable  variations  of  the  true  beam  and  true  target 
parameters.  However,  the  time  constraints  of  doing  this 
study  require  that  only  a  small  variation  of  selected  beam 
parameters  and  target  parameters  be  considered  for 
sensitivity  analysis,  and  similarly  only  a  small  variation 
in  selected  filter  parameters  be  considered  in  robustness 
analysis.  First  this  section  reviews  the  beam,  target  and 
simulation  parameters.  Then  it  presents  the  nominal  values 
for  each  parameter  in  this  study  and  the  sizes  of  variations 
to  be  considered. 

4.3.1  The  Beam  and  Target  Parameters.  This  section 
reviews  the  MMAE  Meer  beam  estimator  and  target  filter 
parameters  used  in  this  study.  The  MMAE  Meer  filter  is 
designed  on  the  basis  of  six  parameters.  They  are  as 
follows : 

rB  is  the  beam  time  constant,  defined  in  Equation  (2-3b). 
This  parameter  is  different  for  each  elemental  Meer 
filter  in  the  bank  structure  of  the  MMAE  algorithm 
g  is  the  square  root  of  the  beam  dynamics  driving  noise 
strength,  g=  (G  Q)  as  defined  in  Equation  (2-3b) 

R  is  the  beam  dispersion  measured  as  the  analog  of  variance 
of  the  Gauss ian-shaped  beam  signal  rate  function  at  the 
detector  surface;  See  Section  (2.2),  Equation  2.7 


An  is  the  noise  arrival  rate  per  length  of  the  detector 
array;  see  Equation  (2-8) 

S  is  the  expected  number  of  signal-induced  events  during  a 
simulation  run;  see  Equation  (2-10) 

D  is  the  elemental  Meer  filter  memory  depth;  see  Section 


L  is  the  length  of  the  one-dimensional  detector  array;  as 
discussed  in  Section  2.2.1/  the  appropriate  minimum  is 
twelve  times  the  beamwidth 

SNR  is  the  signal-to-noise  ratio  as  defined  in  Equation 
(2-10): 


SNR  =  Ifi/(tF-t0) l/I  L) 
=  r  7(  2kR  )  /[XnLJ 


(2-10) 
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Since  fi,  the  number  of  signal -induced  events  expected  during 
one  simulation  run,  (tp-t0),  the  duration  of  one  simulation 
run,  An,  the  noise  arrival  rate  per  length  of  detector,  and 
L,  the  length  of  the  detector,  are  all  independent 
variables,  then  SNR  is  a  dependent  variable.  Therefore,  it 
should  be  determined  from  the  other  parameters.  This  is  in 
contrast  to  previous  thesis  efforts  where  the  noise  rate  was 
made  to  be  a  dependent  variable  using  Equation  (2-10). 

The  three  parameters  associated  with  the  3rd  order 
Gauss-Markov  position  process  target  truth  model  and  target 
Kalman  filter  are  as  follows: 

r<r  is  the  target  time  constant  defined  in  Equation  (3-3) 
is  the  target  dynamics  driving  noise  strength  defined  in 
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Equation  (3-21) 

RT  is  the  target  position  measurement  noise  variance  defined 
in  Equation  (3-8b) 

A  complete  analysis  would  include  the  effects  of  varying 
all  the  parameters  throughout  their  physically  realizable 
ranges.  Time  constraints  prevent  such  a  through  analysis 
and  only  a  few  of  the  parameters  are  considered  in  this 
study.  These  parameters  are  presented  next. 

4.3.2  Nominal  Parameters.  Actual  estimator  and/or 
controller  performance  is  related  to  the  parameters  used  to 
design  it  and  the  conditions  under  which  it  operates.  It  is 
therefore  very  Important  to  know  the  design  parameters  of 
the  filters  and  controllers  as  well  as  the  operating 
parameters  found  in  the  truth  model.  The  following 
parameters  are  provided  as  a  benchmark  for  comparison  with 
later  studies.  The  nominal  values  for  the  beam  and  tracker 
filter  parameters  used  for  the  analyses  are  as  follows: 

Nominal  Beam  Parameters 
7^  =20  sec  D  =  1 

g  =0.2  cm/sec  R  =  0.5  cm^ 

8  =  100  signal  events 

Since  Johnson's  sensitivity  analysis  (the  filter-assumed 
parameter  matches  the  true  parameter  in  sensitivity 
analysis)  found  less  than  1*  variation  in  performance  for 
order  of  magnitude  changes  in  the  beam  time  constant,  the 
nominal  value  for  this  parameter  will  remain  at  20  seconds. 


The  other  beam  parameters  also  remain  the  same  as  those  used 
by  Johnson  so  that  results  will  be  directly  comparable. 

Nominal  Target  Parameters 

Tt  *  10  sec  R  =  0.5  cm^ 

2  5 

Q  #1  =  0.01  cm  /sec 

Q  #2  =  0.07  cm2/sec5 

Q  13  =  1.00  cm2/sec^ 

The  three  values  of  QT  correspond  to  the  three  values  chosen 
by  Johnson  [9:691  for  adaptively  modelling  the  target  noise 
strength.  Johnson  describes  these  three  values  as  modelling 
the  trajectories  from  the  straight  and  level  flight,  to  the 
maximum  manned-vehicle  g-limit  of  10  g's. 

The  simulation  Parameters 

Zicker  studied  the  number  of  Monte  Carlo  runs  necessary 
for  the  sample  statistics  to  converge  to  a  consistent  value 
(321.  Based  on  Zicker's  work,  the  number  of  Monte  Carlo 
runs,  N,  is  set  at  200.  This  large  number  is  necessary  due 
to  the  sparse  number  of  photon  events. 

Because  the  fastest  transients  at  the  nominal  condition 
have  a  time  constant  of  ten  seconds  (rather  benign 
dynamics),  the  Kalman  filter  and  the  feedback  control  sample 
period  will  be  set  to  one  second.  This  easily  satisfies  the 
Shannon  sampling  theorem  which  states  that  the  sampling  rate 
should  be  at  least  twice  the  highest  signal  frequency 
content  of  interest  (33:87-881.  One  second  also  provides 
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the  extra  margin  suggested  by  the  engineering  rule  of  thumb, 
of  sampling  five  to  ten  times  faster  than  the  highest 
frequency  used  [33:4971. 

As  previously  noted,  the  detector  length  is  constrained  by 

the  beamwidth.  However,  the  choice  of  detector  length  is 

further  constrained  by  the  scaling  chosen  and  computer 

limitations  for  exponent  arguments.  The  value  of  20  cm  for 

2 

L  allows  adequate  size  for  beam  dispersions  from  .3  cm  to 
2 

.7  cm  ,  while  maintaining  exponent  arguments  within  computer 
defined  limits.  This  L  is  adequate  for  this  study  and  it  is 
desirable  to  maintain  the  same  scale  with  past  studies, 
however  future  studies  may  wish  to  apply  scaling  factors  to 
the  problem  to  avoid  limitations  on  exponent  argument 
limits . 

To  provide  a  Xn  that  is  consistent  with  the  values  used 
in  the  past  (  A.n  corresponding  to  SNR  -  20  and  L  =  10),  when 
the  mean  of  signal  rate  was  1  event/sec,  An  is  set  to  0.005 
events/sec/cm.  Since  the  detector  length  is  now  20  cm,  the 
equivalent  SNR  is  now  10.  Having  a  constant  noise  rate  per 
unit  length,  rather  than  one  dependent  on  the  SNR,  L, 
expected  number  of  signal  events,  and  duration  of  the  run, 
will  provide  more  realistic  results  for  changes  in  signal 
event  rate,  and  detector  length. 

4.3.3  Analysis  Parameters.  Since  Johnson  experienced 


significant  degradation  in  controller  tracking  performance 
with  only  2\  variation  in  true  beam  time  constant  from  the 


filter  assumed  value  (9:1181/  this  Is  the  focal  point  for 
robustness  analysis  of  the  controller.  The  2%  variation  is 
also  the  starting  point  for  the  MMAE  estimator  performance 
evaluation.  Larger  variations  of  the  true  beam  time 
constant  from  the  nominal,  filter-assumed  value  may  be 
studied  if  beam  estimator  performance  or  controller  tracking 
performance  at  2%  variation  warrants  further  study.  Johnson 
found  no  sensitivity  problems  with  variations  of  the  beam 
time  constant  nor  significant  sensitivity  or  robustness 
problems  with  other  beam  parameters  (91.  The  inclusion  of 
the  MMAE  structure  within  the  controller  may  result  in  a 
change  to  the  sensitivity  or  robustness  of  the  closed  loop 
system.  Therefore,  sensitivity  and  robustness  analyses  are 
performed  for  variations  in  the  beam  dynamics  driving  noise 
strength  and  beam  dispersion.  The  variation  considered  will 

be  one  order  of  magnitude  above  and  below,  with  a  nominal 
2 

value  of  0.2  cm  /sec  for  the  beam  dynamics  driving  noise 

2 

strength  and  a  nominal  value  of  0.5  cm  for  the  beam 
dispersion.  This  will  allow  comparisons  to  be  made  with 
Johnson's  results  (91.  The  next  chapter  will  describe  the 
results  obtained  for  simulation  runs  based  on  the  nominal 
parameters  given  in  this  section. 


4 • 4  Summary 

The  purpose  of  this  chapter  has  been  to  explain  the 
tools  and  methods  used  to  evaluate  estimator  performance  of 
the  MMAE  Meer  filter  and  to  evaluate  the  controller 
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performance  of  the  MMAC  based  on  weighting  matrices 
generated  through  implicit  model  following  techniques.  The 
Monte  Carlo  simulation  provides  the  most  viable  method  of 
evaluating  the  estimator's  or  controller's  performance. 

SO FE  and  SOFEPL  provide  the  basic  method  of  computing  the 
filter  state  estimation  error  statistics  and  controller 
tracking  error  statistics,  which  allow  us  to  evaluate 
controller  designs.  Four  of  the  performance  analyses: 
controller  performance  with  MMAE  Meer  filter,  controller 
performance  with  reduced  depth  MMAE  Meer  filter,  controller 
performance  with  the  expanded  bank  MMAE  Meer  filter,  and  the 
time-varying  parameter  analyses,  determine  the  capabiltiies 
of  the  MMAE  Meer  filter  to  handle  off-nominal  beam  time 
constants  in  a  closed-loop  system.  Three  of  the  performance 
analyses:  MMAC  performance  with  weighting  matrices  generated 
from  implicit  model  following  techniques,  MMAC  sensitivity 
investigation,  and  MMAC  robustness  analysis,  evaluate  the 
capabity  of  weighting  matrices  based  on  implicit  model 
following  techniques  to  provide  additional  robustness  to  the 
closed-loop  system  incorporating  the  MMAE  Meer  filter  while 
maintaining  acceptable  performance  at  design  conditions. 

The  results  of  the  seven  analyses  are  in  the  next  chapter. 


V.  RESULTS  AND  ANALYSIS 

The  results  of  the  eight  Monte  Carlo  analyses  discussed 
in  Chapter  4  are  presented  in  this  chapter.  The  first  five 
analyses  are  on  the  MMAE  Meer  filter  operating  in  open  loop 
as  a  beam  estimator.  The  last  three  analyses  are  on  the 
MMAE  Meer  filter  operating  in  closed  loop.  In  this  mode, 
the  MMAE  Meer  filter  provides  state  estimates  to  a 
controller  which  uses  the  beam  estimate,  combined  with 
target  state  estimates,  to  generate  a  control  input  to  point 
the  beam  at  the  target.  The  controller  used  is  the  three- 
element  MMAC  developed  by  Johnson  [9]  with  constant  gains 
based  on  the  nominal  beam  time  constant. 

The  first  analysis  demonstrates  the  MMAE  Meer  filter  can 
provide  precise  estimates  of  the  beam  location  when  the  beam 
time  constant  is  uncertain.  The  second  analysis  shows  that 
the  same  MMAE  Meer  Filter  does  not  provide  a  precise 
adaptive  estimate  of  the  beam  time  constant  when  the  beam  is 
not  in  a  closed  loop.  The  third  analysis  indicates  use  of 
an  elemental  Meer  filter  depth  of  one  is  adequate  for  state 
estimation  purposes.  The  fourth  MMAE  Meer  filter  analysis 
shows  that  the  filter  can  perform  well  in  the  presence  of  a 
linearly  changing  or  a  sinusoidally  varying  true  beam  time 
constant.  The  final  MMAE  Meer  filter  analysis  indicates 
that  a  simple  bank  of  3  Meer  filters  performs  well  for  a 
broad  range  beam  time  constants. 


Operation  of  the  MMAE  Meer  filter  in  a  closed  loop 
tracker  results  in  a  change  of  performance  for  the  MMAE  Meer 
filter  from  that  seen  in  the  open  loop  analysis.  To  ensure 
that  signal  events  occur  on  the  array  while  the  beam  is 
tracking  the  target,  two  methods  of  keeping  the  beam  over 
the  detector  were  implemented  and  analyzed.  The  first  is, 
conceptually,  the  rotation  of  the  NPB  toward  the  target. 

The  second  is,  conceptually,  the  sweeping  of  the  detector 
underneath  the  beam.  Subsequent  analyses  employ  the  second 
approach.  The  performance  of  the  MMAE  Meer  filter  is  then 
reanalyzed  for  closed  loop  operation.  Along  with  this 
analysis,  tracking  performance  baselines  are  provided 
against  which  subsequent  controller  designs  may  be  compared. 

Time  constraints  prevented  the  analyzing  the  performance 
of  the  nine-element  MMAC  based  on  implicit  model  following 
weighting  matrices.  However,  performance  of  this  MMAC  is 
presented,  using  the  state  and  control  weightings  of: 

xn  =  100,  x55  =  10,  U  =  1 

where  these  are  the  quadratic  cost  weighting  matrices  on  the 
target  tracking  error,  the  pseudointegral  state  and  control, 
respectively.  This  allows  a  comparison  of  this  MMAC  with 
the  baselines  and  past  controller  designs. 

5.1  MMAE  Meer  Filter  State  Estimation  Performance 


The  first  analysis  indicates  the  MMAE  Meer  filter  can 


provide  precise  estimates  of  the  beam  location  when  the  beam 
time  constant  is  uncertain.  The  study  uses  an  MMAE  Meer 
filter  composed  of  three  elemental  Meer  filters.  The  beam 
time  constants  for  Meer  filters  1,  2,  and  3  are  20.0  20.8 
and  19.2  respectively,  representing  a  nominal  beam  time 
constant  plus  and  minus  four  percent.  Each  Meer  filter  is 
provided  the  true  value  of  the  beam  dispersion  and  the 
dynamics  driving  noise  strength.  Since  all  the  Meer  filter 
and  MMAE  Meer  filter  variables  have  to  be  sampled  at  regular 
intervals  for  ensemble  averaging  of  the  Monte  Carlo  run,  it 
is  necessary  to  propagate  the  state  estimates  of  each 
individual  Meer  filter  and  the  MMAE  MMAE  state  estimate. 
Thus,  it  is  necessary  to  select  a  beam  time  constant  to 
propagate  the  MMAE  Meer  filter  state  estimate  to  each  sample 
time.  The  choices  considered  were  the  adaptive  beam  time 
constant  of  the  MMAE  filter  and  the  nominal  time  constants 
of  the  three  filters.  An  initial  run,  using  the  adaptive 
beam  time  constant  to  propagate  the  adaptive  state  estimate 
to  the  sample  times,  was  made  with  true  beam  time  constant 
of  20.4  seconds.  This  run  indicated  that  the  estimate  of 
the  beam  time  constant  was  poor.  Therefore,  the  nominal 
beam  time  constant  is  used  to  propagate  the  state  estimate 
in  subsequent  runs  of  the  MMAE  Meer  filter  in  open  loop.  In 
addition,  the  Meer  filter  initial  state  variances  were 
adjusted  from  the  past  value  of  .1  to  .4.  The  latter  value 
is  the  mean  square  value  for  centroid  oscillations, 
representing  the  steady  state  motion  of  the  beam  centroid. 


as  derived  from  Equation  (2-16): 

PSS  ■  1/2  tb  92  <*-!> 

The  rms  errors  for  various  values  of  the  true  beam  time 
constant  are  in  Table  5.1.  As  the  table  indicates,  the  RMS 
error  for  lower  beam  time  constants  is  smaller.  This  trend 
is  due  to  the  lower  steady  state  oscillation  of  the  beam  as 
given  by  Equation  (5-1).  Figure  5.1  is  the  MMAE  Meer  mean 
error  plus  and  minus  one  standard  deviation  time  history  for 
true  beam  time  constant  of  20.4  seconds.  However,  Figure 
5.1  is  also  representative  of  the  elemental  Meer  filter 
performance,  regardless  of  the  true  value  of  the  beam  time 
constant.  The  central  jagged  line  is  the  mean  error  history 
of  the  MMAE  estimate  from  the  true  beam  position.  The 
smooth  dashed  lines  are  the  envelope  of  plus  and  minus  the 
filter's  predicted  error  standard  deviation  or  square  root 
of  the  filter-computed  variance.  The  outer  jagged  lines  are 
the  mean  error  ±  one  standard  deviation  envelope.  What  is 


Table  5.1  MMAE  Meer  State  Estimate  Performance 


2 

R  =  .  5  cm  g  = 

1/2 

.2  cm/sec  Depth  =  3 

True  beam 

RMS  error 

time  constant 
(sec) 

(cm) 

20.4 

.377023 

20.0 

.376193 

19.6 

.374878 

significant  in  the  plot  is  that  the  mean  error  is  near  zero 
and  the  filter  predicted-error  standard  deviation  remains 
near  the  true  standard  deviation  of  the  error.  The 
performance  of  the  individual  Meer  filters  is  given  in  Table 
5.2.  As  noted,  there  is  very  little  variation  in  individual 
Meer  filter  performance  for  changes  in  the  beam  time 
constant.  Thus,  the  Meer  filter  has  little  difficultly  in 
estimating  the  beam  location  with  a  mismodelled  beam  time 
constant  for  an  uncontrolled  beam.  As  previously  discussed, 
the  RMS  errors  drop  as  the  true  beam  time  constant  is 
reduced  due  the  lower  steady  state  oscillations  of  the  beam. 
Also  note  that  the  Meer  filter  with  the  highest 
filter-assumed  beam  time  constant  has  the  lowest  RMS  error 
of  the  three  Meer  filters.  Again  this  is  related  to  the 
steady  state  error  variance.  The  larger  filter-assumed  beam 
time  constant  results  in  a  larger  filter  computed  variance. 
This  larger  variance  results  in  a  slightly  larger  gain 
placed  on  incoming  measurements.  This  results  in  a  better 
estimate  of  the  beam  location  and  consequently  lower  RMS 
error.  Table  5.3  is  presented  to  show  the  similarity  in 
performance  between  the  filters.  No  other  conclusion  should 
be  drawn  from  Table  5.3  as  residual  performance  for  the  Meer 
filter  involves  both  good  tracking  of  the  beam  (resulting  in 
lower  residuals)  and  good  rejection  of  noise  events 
(resulting  in  high  residuals). 

Figure  5.2  shows  the  residual  mean  history  of  an 
elemental  Meer  filter.  As  before,  this  figure  is 


Table  5.2  Individual  Meer  State  Estimate  Performance 
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TB=  20.0 

tb=20 . 8 
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.377067 

.376949 

.377254 
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20.0 

.376225 

.376138 

.376378 

19.6 

.374901 

.374893 

.374969 

Table  5.3  Meer  Filter  Pseudo-Residual  Performance 


2 

R  =  . 5  cm  g  = 

1/2 

.2  cm/sec  Depth 

=  3 

True  beam 

RMS 

residual 

time  constant 

(cm) 

(sec) 

Adaptive 

Meerll 

Meer#2 

Meerf 3 

Tg®  20.0 

Tb=20.8 

7g=  19.2 

20.4  .129923 

.129914 

.129956 

.129869 

20.0  .129907 

.129900 

.129942 

.129854 

19.6  .129921 

.129913 

.129955 

.129869 

representative  of  the  residual  history  of  any  of  the  Meer 
filters.  Of  importance  here  is  that  the  mean  of  the 
residuals  remains  within  the  filter  computed  tier  bounds, 
shown  on  the  plot  as  dashed  lines.  Of  equal  importance  is 
the  fact  that  the  results  vary  so  little  with  changes  in 
beam  time  constant.  It  should  be  noted  that  these  are 
pseudo-residuals.  As  with  the  state  estimates,  it  is 
necessary  to  collect  the  residuals  at  a  uniform  sampling 
interval  to  allow  ensemble  averaging.  The  sampling  Interval 
is  one  second  to  match  the  average  number  of  expected  signal 
events.  However,  Meer  filter  measurements  occur  at  random 
intervals.  If  two  events  occur  during  a  sample  period,  then 
only  the  residual  from  the  second  event  is  processed.  If  no 
events  occur  during  the  sample  period,  then  the  previous 
residual  is  processed  again.  Thus,  some  residuals  are  not 
collected  for  processing  while  a  few  may  be  processed  more 
than  once.  Therefore,  these  pseudo-residuals  are  only  used 
to  show  gross  trends  in  the  data.  As  the  performance  is 
remarkably  good,  no  additional  Meer  filters  were 
incorporated  into  the  MMAE  structure.  The  next  results  are 
on  the  adaptive  beam  time  constant. 

5.2  MMAE  Meer  Filter  Parameter  Estimation  Performance 


conducted  In  conjunction  with  the  first  analysis.  The  beam 
time  constants  for  each  Meer  filters  are  again,  20.0,  20.8 
and  19.2  seconds.  The  average  values,  over  the  100  second 
interval  of  the  runs,  for  the  mean  value  of  the  adaptive 
beam  time  constant  and  the  average  value  of  the  filter 
computed  variance  of  the  adaptive  estimate  are  in  Table  5.4. 
Of  importance  here  is  that  once  again  there  is  little 
variation  in  the  performance  for  changes  in  the  beam  time 
constant.  The  difference  in  the  variance  for  the  run  with 
true  beam  time  constant  of  20.4  seconds  is  due  to  an  error 
in  setting  the  initial  value  of  the  variance  to  0  for  that 
run.  As  one  might  expect  from  the  results  of  the  first 
study,  the  parameter  estimation  is  very  poor.  This  is  due 
to  the  small  difference  in  the  performance  between  the  Meer 
filters.  Under  these  conditions  one  would  expect  the 
adaptive  parameter  estimate  will  tend  toward  the  filter  with 
the  smallest  |A^|,  where  A  ^  is  defined  in  Section  2.4 
Equation  (2-33)  [17:1331.  Therefore,  one  would  expect  the 

Table  5.4  MMAE  Meer  Parameter  Estimate  Performance 

2  1/2 

R  =  . 5  cm  g  =  . 2  cm/sec  Depth  *  3 

True  beam  Adaptive  beam  Filter  computed 

time  constant  time  constant  Variance 

(sec)  (sec) 

20.4  20.1022  .408476 

20.0  20.1021  .432341 

19.6  20.1019  .432323 


estimate  would  tend  toward  the  parameter  estimate  associated 
with  the  filter  with  the  smallest  beam  time  constant  as  that 
would  result  in  the  lowest  steady  state  variance  as  given  by 
Equation  (5-1).  However,  for  the  Meer  filter,  this  is 
offset  by  the  improved  residual  performance  of  the  filter 
with  the  largest  beam  time  constant.  The  effect  of  this 
slight  improvement  in  residual  performance  is  higher 
probabilities  for  the  filter  with  the  largest  l*kl‘  Figure 
5.3  depicts  the  time  history  of  the  adaptive  beam  time 
constant.  Figures  5.4  through  5.6  show  the  individual 
filter  probabilities  for  a  beam  time  constant  of  20.8 
seconds.  Again  these  figures  are  representive  of  the 
performance  with  any  selection  of  true  beam  time  constant. 

As  with  the  pseudo-residuals,  the  adaptive  beam  time 
constant  and  the  probabilities  are  calculated  at  measurement 
times,  but  are  recorded  at  one  second  intervals  for  ensemble 
averaging  of  the  Monte  Carlo  data.  Figure  5.3  might  at 
first  appear  to  indicate  that  the  beam  time  constant  is 
slowly  approaching  the  correct  value.  Actually,  it  is 
slowly  moving  toward  the  value  of  20.8  seconds  of  Meer 
filter  82,  as  was  explained.  To  demonstrate  this,  a  Monte 
Carlo  run  extended  to  a  final  time  of  300  seconds  is 
displayed  in  Figure  5.7.  The  adaptive  beam  time  constant  is 


prevented  from  actually  reaching  20.8  seconds  due  to  a  lower 
bound  on  probability  limit  on  each  estimator  of  .1.  Thus, 
the  highest  adaptive  estimate  possible  with  a  lower 
probability  bound  of  .1  is  only  20.56.  The  0.1  bound  was 
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selected  to  match  the  lower  bound  set  by  Johnson  [9]  for  the 
Kalman  filter  probabilities  in  the  MMAC. 

The  intent  of  the  .1  probability  bound  is  to  prevent 
delays  in  switching  between  a  controller  using  a  Kalman 
filter  based  on  a  high  Q  value  for  a  highly  dynamic  target 
and  a  controller  based  on  a  Kalman  filter  with  a  low  Q  value 
for  a  benign  target.  However ,  a  bound  of  .1  is  generally 
considered  to  be  excessive  even  for  this  purpose.  Later 
runs  use  a  revised  Meer  filter  probability  bound  of  0.0005 
as  a  more  appropriate  value  for  rather  benign  beam  dynamics. 
This  value  allows  the  adaptive  value  to  approach  the  upper 
and  lower  limits,  without  allowing  a  filter's  probability  to 
drop  to  zero.  If  the  filter's  probability  ever  went  to 
zero,  then  it  would  no  longer  have  an  impact  on  the  adaptive 
state  or  adaptive  parameter  estimate,  regardless  of  the 
actual  beam  time  constant.  A  lower  bound  that  is  too  small 
might  delay  the  response  to  a  dynamic  beam  time  constant. 

For  the  beam  in  open  loop,  the  lower  probability  bound  has 
little  impact  for  the  short  term  case,  as  there  is  little 
variation  in  the  Meer  filter  probabilities  as  indicated  by 
Figures  5.4,  5.5  and  5.6.  Later  runs  with  the  new  lower 
bound  verified  this  conclusion. 

5.3  Simplified  MMAE  Meer  Filter  State  Estimation  Performance 

The  third  analysis  indicates  use  of  an  elemental  Meer 
filter  depth  of  one  is  adequate  for  state  estimation 
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purposes.  In  this  study,  special  coding  is  used  to 
Implement  the  simplified  update  routine  described  in  Section 
2. 5. 2. 2  Equation  (2-46).  In  addition,  the  Meer  filter 
propagating  routine  is  streamlined  for  the  Depth  =  1  case. 

As  indicated  in  Table  5.5,  the  MMAE  Meer  state  estimation 
performance  showed  little  difference  at  a  depth  =  1  compared 
to  depth  =  3.  This  is  not  surprising  when  the  effect  of  the 
filter  depth  on  the  Merge  method  of  Section  2. 5. 2. 2  is 
considered.  The  result  of  increasing  the  filter  depth  Is  to 
extend  the  amount  of  time  that  a  measurement  affects  the 
Meer  filter  estimate  (and  the  MMAE  Meer  estimate).  As  such, 
noise  measurements  continue  to  have  an  effect  on  the 
estimate  until  eliminated  from  the  filter.  The  purpose  of 
the  Merge  method  is  to  retain  a  greater  portion  of  the 
measurement  history  in  collapsing  the  filter  hypothesis  tree 
than  can  be  achieved  by  acting  on  a  decision  to  discard  half 
of  the  measurement  hypotheses  as  in  the  "Best  Half"  method 
of  Section  2. 5. 2.1.  So  a  measurement  retains  an  indirect 


Table  5.5  State  Estimate  Performance  and  Computional 
Loading  For  Depths  of  1  and  3 


2  cm/sec 
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Depth  =  3 
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Impact  on  the  state  estimation  after  leaving  the  filter. 
Thus,  a  depth  of  one  still  retains  much  of  the  measurement 
history  while  also  minimizing  the  impact  of  noise  events. 

As  shown  in  the  table,  there  is  a  great  deal  of  time 
variation  between  runs.  This  variation  is  due  to  computer 
loading  at  the  time  of  the  run.  Most  runs  were  conducted 
during  evening  hours  when  computer  loading  is  light.  A 
notable  exception  is  the  run  with  true  beam  time  constant  of 
19.6  and  Depth  of  1,  that  was  run  mid-morning.  Comparison 
of  runs  made  at  similar  times,  indicate  there  may  be  some 
slight  reduction  in  computer  time  when  using  a  depth  of  1. 
Two  reasons  account  for  such  a  small  change  in  computation 
time  for  a  change  in  depth.  First,  the  code  was  altered  to 
provide  greater  sensitivity  to  all  depths.  As  an  example, 
the  probability  calculation  associated  with  Equation  (2-40) 
was  performed  for  256  elemental  Snyder -Fishman  filters 
regardless  of  depth  in  the  software  of  past  studies.  The 
use  of  3  Meer  filters  would  have  tripled  this  number, 
however,  the  code  for  this  study  was  altered  so  only  24 
probability  calculations  were  required  at  a  depth  of  3  for 
all  3  Meer  filters.  This  is  only  8  times  the  number  of 
probability  calculations  necessary  for  3  Meer  filters  at  a 
depth  of  3.  Second,  as  indicated  in  the  SOFE  manual 
(24:831,  the  routines  involved  in  propagating  the  truth  and 
Kalman  filter  states  via  the  Kutta-Merson  integration 
routine  are  called  most  often.  This  is  a  substantial 
portion  of  the  computer  load.  The  elemental  Snyder-Fishman 


filter  propagation  is  accomplished  separately,  using  a  one 
dimensional  state  transition  matrix.  Changes  in  the  Meer 
filter  depth  will  not  have  a  large  impact  unless  the  signal 
rate  is  high.  Even  so,  the  computer  loading  will  still  be 
dominated  by  the  integration  routine  as  it  propagates  the 
truth  and  filter  states  to  the  Poisson  event  times,  as  well 
as,  noise  injection  and  Kalman  update  periods.  To  calculate 
the  loading  more  specifically  due  to  the  Meer  filter 
operations  one  could  count  the  number  of  arithmetic 
operations  per  Meer  cycle  for  changes  in  depth.  However,  as 
the  Meer  filter  is  one  dimensional  (thus  one  state  and  one 
variance),  the  number  of  arithmetic  operations  will  increase 
by  2  raised  to  the  Depth  (D)  power.  Thus,  the  Meer  filter 
computation  load  for  a  depth  of  3  as  compared  to  a  depth  of 
1  is  8  times  greater.  If  the  Kalman  filter  were  also 
propagated  by  a  state  transition  matrix,  the  same  approach 
could  be  used  to  establish  the  relative  loading  between  the 
Kalman  and  Meer  filters.  However,  that  analysis  was  not 
performed  in  this  study.  Based  on  the  results  in  Table  5.5, 
the  remaining  analyses  were  conducted  at  Depth  =  1,  with  the 
Simplified  MMAE  Meer  filter  code. 


5 . 4  MMAE  Meer  Performance  With  Variable  Beam  Time  Constant 


This  study  consisted  of  altering  the  beam  time  constant 
of  the  truth  model  without  Informing  the  MMAE  Meer  filter. 
As  before  the  beam  time  constants  of  filters  1,  2,  and  3 
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were  20.0,  20.8,  and  19.2  seconds  respectively.  The  true 
beam  time  constant  was  varied  from  an  initial  value  of  19.2 
seconds  by  adding  a  factor  of  0.016t  to  the  beam  time 
constant  in  the  propagation  routine.  These  values  were 
chosen  to  span  the  range  of  beam  time  constants.  The  result 
is  a  ramping  beam  time  constant.  This  was  then  repeated  for 
a  descending  ramp,  with  an  initial  true  beam  time  constant 
of  20.8  and  slope  of  -0.016t.  Table  5.6  has  the  state 
estimation  performance  under  these  conditions.  Figure  5.8 
depicts  the  mean  adaptive  beam  time  constant  for  the  case  of 
the  true  beam  time  constant  descending  linearly  from  a  value 
of  20.8  to  19.2.  As  indicated  in  Figure  5.8,  the  adaptive 
beam  time  constant  is  also  ramping.  However,  the  positive 
slope  matches  that  seen  in  Figure  5.3,  rather  than  the 
negative  input  slope  of  -0.016t.  Again,  this  figure  is 
characteristic  of  the  results  obtained  for  the  beam  time 
constant  with  a  positive  ramp  and  for  later  results  with  a 

Table  5.6  State  Estimate  Performance 
With  Variable  Beam  Time  Constant 


2  1/2 

R  =  . 5  cm  g  *  .2  cm/sec  Depth  =  1 


True  beam 

RMS  error 

time  constant 

(cm) 

(sec) 
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19. 2+.016t 

.375747 

20.8- . 016t 

.375646 

20 . 0+ . 8sin(omega  t) 

.375633 

omega  =  2n/100 


sinusoidally  varying  beam  time  constant. 


Next,  the  beam  time  constant  was  varied  in  a  sinusoidal 
manner  by  adding  a  factor  of  0.8  sin(omega  t).  The  initial 
beam  time  constant  was  again  20.0  seconds.  A  magnitude  of 
0.8  was  chosen  to  match  the  range  of  beam  time  constants  of 
the  filters.  The  argument  of  the  sin  function,  omega,  %/as 
chosen  to  be  2x/100  to  provide  a  complete  cycle  in  100 
seconds.  The  results  are  also  in  Table  5.6  and  again  are 
similar  to  previous  results. 

5.5  MMAE  Meer  Filter  Range  With  Three  Meer  Filters 

The  final  MMAE  Meer  filter  analysis  indicates  that  a 
simple  bank  of  3  Meer  filters  performs  well  for  a  broad 
range  of  beam  time  constants.  As  the  MMAE  Meer  performance 
in  previous  analyses  seemed  insensitive  to  variations  in  the 
beam  time  constant,  the  author  chose  filter  1,  2,  and  3  beam 
time  constants  of  20,  200,  and  2  seconds,  respectively,  for 
this  test.  This  choice  not  only  represents  two  orders  of 
magnitude  spread  in  the  beam  time  constant,  it  stresses  the 
sampling  rate.  The  true  beam  time  constant  was  set  at  a 
constant  value  of  70  seconds,  as  this  is  about  halfway 
between  20  and  200  on  a  log  scale.  Figure  5.9  is  the  mean 
error  history  +/-  1  standard  deviation  for  the  adaptive 
state  estimate.  This  figure  is  also  representative  of  the 
performance  of  the  Meer  filters  with  beam  time  constants  of 
20  seconds  and  200  seconds.  Figure  5.10  is  the  performance 


of  the  Meer  filter  with  beam  time  constant  of  2  seconds.  Of 
significance  in  Figure  5.9  is  that  the  mean  error  is  near 
zero  and  the  filter-computed  standard  deviation  is  close  to 
the  actual  standard  deviation.  In  Figure  5.10  we  see  that 
the  mean  error  does  not  remain  near  zero  and  the  filter 
predicted  standard  deviation  is  much  less  than  the  actual 
standard  deviation.  This  is  an  indication  of  degraded 
filter  performance  due  to  a  poor  assumed  beam  time  constant. 
Note  that  in  Table  5.7,  the  RMS  error  for  the  Meer  filter 
with  beam  time  constant  of  200  is  lower  than  the  RMS  error 
for  the  other  filters.  This  is  due  to  the  gain  placed  on 
the  measurements  through  Equation  (2-22)  that  is  part  of  the 
weighting  applied  to  the  residuals  in  Equation  (2-46).  The 
influence  of  the  large  beam  time  constant  on  Equation  (2-22) 
is  inferred  through  Equation  (5-1):  see  the  last  column  of 
Table  5.7.  Also  note  that  the  RMS  error  is  slightly  lower 
than  the  filter  predicted  standard  deviation.  Thus,  this 


Table  5.7  Stressed  Filter  Performance 
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filter  Is  performing  well.  In  contrast,  the  filter 
predicted  standard  deviation  of  the  Meer  filter  with  beam 
time  constant  of  2  is  much  lower  than  the  RMS  error.  This 
filter  is  operating  improperly.  For  this  filter,  the 
variance  in  steady  state  is  much  smaller  as  inferred  by 
Equation  (5-1).  Thus,  measurements  in  a  filter  with  a  low 
beam  time  constant  receive  correspondingly  lower  weight 
through  Equation  (2-22)  and  the  filter  relies  on  its 
internal  dynamics  model.  When  the  filter  dynamics  model  is 
incorrect  (as  with  this  filter)  the  state  estimation 
suffers.  The  degraded  state  estimate  results  in  higher 
residuals  for  actual  signal  events.  Higher  residuals  result 
in  lower  weighting  on  the  measurements  through  Equations 
(2-7),  (2-38),  (2-39)  and  (2-46).  Consequently,  the  filter 
that  underestimates  the  beam  time  constant  diverges.  Note 
that  for  the  divergent  filter  (Meer  filter  #3),  the 
predicted  value  is  approximately  equal  to  the  steady  state 
variance  defined  in  Equation  (5-1)  for  a  filter  based  on  a 
beam  time  constant  of  2  seconds.  This  indicates  that  the 
measurements  are  being  discarded  as  noise.  As  the  higher 
residuals  of  the  divergent  filter  also  result  in  a  low 
probability  weight  on  the  state  estimates  from  that  filter, 
the  MMAE  Meer  is  able  to  reject  the  influence  of  a  filter 
that  underestimates  the  beam  time  constant.  Although  the 
filter  computed  variance  of  the  Meer  filter  with  beam  time 
constant  of  20  seconds  is  lower  than  the  steady  state 
variance  Indicated  by  Equation  (5-1)  for  a  beam  time 


constant  of  20  seconds,  its  square  root  is  lower  than  the 
RMS  error.  Therefore,  the  filter  is  making  use  of  the 
incoming  measurements;  however,  the  weighting  is  not 
sufficient  to  track  the  beam  precisely.  The  RMS  error  of 
the  adaptive  state  estimate  also  exceeds  its  filter-computed 
error  standard  deviation.  This  is  not  surprising  as  the 
adaptive  state  estimate  for  this  run  uses  the  nominal  beam 
time  constant  of  20  seconds  rather  than  the  true  beam  time 
constant  of  70  seconds  to  propagate  to  the  measurement  time. 
Also,  this  run  was  conducted  prior  to  the  setting  of  the 
Meer  filter  lower  probability  bound  from  .1  to  0.0005. 
Therefore,  the  divergent  filter  impacts  the  adaptive  state 
estimate  to  a  greater  degree  than  would  normally  be 
expected . 

The  exceptional  performance  of  the  Meer  filter  with  beam 
time  constant  of  200  becomes  a  problem  for  the  adaptive 
estimation  process.  The  adaptive  parameter  estimate  is 
quickly  pulled  past  the  correct  beam  time  constant  of  70 
seconds  toward  the  filter  with  beam  time  constant  of  200,  as 
shown  in  Figure  5.11.  The  upper  limit  on  the  adaptive 
estimate  due  to  the  0.1  lower  probability  bound  is  162 
seconds.  However,  the  influence  of  the  Meer  filter  with 
beam  time  constant  of  20  seconds  results  in  a  mean  adaptive 
beam  time  constant  of  125  seconds  (averaged  over  the  final 
50  seconds  of  the  simulation).  This  is  a  poor  adaptive 
estimate,  but  not  a  surprising  result  considering  the  slow 
trend  toward  the  Meer  filter  with  the  highest  time  constant 
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in  previous  analyses.  Thus,  Meer  filters  with  large  beam 
time  constants  may  dominate  the  adaptive  beam  time  constant 
estimation  process.  It  may  be  possible  to  reduce  this 
influence  by  decreasing  the  dynamics  driving  noise  term  g  as 
used  in  Equation  (5-1).  However,  such  detuning  may  have 
undesirable  impacts  on  state  estimation. 


5 . 6  MMAE  Meer  Filter  In  Controller  Looi 


The  results  of  the  previous  sections  indicate  that  the 
Meer  filter  itself  is  rather  robust  with  respect  to  the  beam 
time  constant.  Then  one  asks  why  did  Johnson  [9]  experience 
instability  in  the  Meer  filter  when  operated  in  the  closed 
loop?  This  section  explores  that  question  and  follows  the 
attempts  to  stabilize  the  MMAE  Meer  filter  in  closed  loop 
operation . 


5-6-1  Beam  Following  Detector.  The  first  step  in 
exploring  the  closed  loop  operation  was  another  review  of 
the  code  to  ensure  there  were  no  errors  in  the  control  loop. 
This  review  disclosed  a  significant  error.  In  the 
Implementation  of  the  third  order  Gauss-Markov  position 
process,  no  provision  was  made  to  ensure  that  the  beam  would 
remain  over  the  detector.  As  the  position  of  the  target 
went  to  infinity,  the  controller  forced  the  beam  off  the 
array  to  keep  it  on  the  target.  Meer  indicated  [21:1811, 
that  the  probability  of  a  noise  event  outside  the  detector 
is  zero.  Therefore,  the  program  set  the  probability  of  an 
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assumed  noise  event  to  zero  and  the  signal  event  probability 
became  one,  whenever  an  event  occurred  off  the  array.  Meer 
did  take  into  consideration  the  finite  size  of  the  detector 
on  the  generation  of  signal  events  [21:149].  However,  in 
Meer's  problem  it  was  appropriate  to  set  the  array  size  much 
larger  than  the  beamwidth  and  assume  the  probability  of  a 
signal  off  the  array  as  negligible.  For  pointing  of  the 
beam  at  a  highly  dynamic  target,  such  assumptions  have  to  be 
reconsidered.  As  the  initial  target  position  and  velocity 
were  set  to  0  cm  and  5  cm/sec,  the  10  cm  detector  used  in 
previous  studies  was  generally  left  behind  in  the  first  two 
seconds.  This  effectively  reduced  the  Meer  filter  to  a 
Snyder-Fishman  filter  with  respect  to  signal  events.  The 
author  believes  this  error  cannot  truly  account  for  the 
Instability  of  the  estimator  in  the  closed  loop  exhibited  in 
Johnson's  research  [9],  but  it  may  have  aggravated  it.  As  a 
Snyder-Fishman  filter,  the  Meer  filter  would  weight  each 
incoming  signal  event  as  an  actual  event.  There  would  be  no 
deweighting  of  the  event  due  to  a  large  residual.  However, 
a  noise  event  would  still  occur  on  the  array  so  that  the 
Meer  filter  would  deweight  the  residual  based  on  the  size  of 
the  residual  through  Equations  (2-7),  (2-39)  and  (2-46). 
Therefore,  noise  events  had  neglible  impact,  while  the 
filter  was  given  full  benefit  of  the  information  in  signal 
events.  Basically,  the  Meer  filter  was  being  told  which 
events  were  noise  induced  and  which  events  were  signal 
induced.  Therefore,  this  error  cannot  account  for  the 


Instability  problem.  However,  it  probably  did  have  Impacts 
on  the  signal  to  noise  (SNR)  sensitivity  and  robustness 
studies  done  by  Johnson  [91  and  his  sensitivity  studies 
involving  variations  in  the  number  of  expected  signal  events 
over  the  run. 

The  author  considered  three  approaches  to  maintaining 
the  beam  on  the  array.  The  first  approach  was  to  increase 
the  slzf  of  the  detector.  However,  this  would  require  a 
length  of  detector  of  several  meters  and  therefore  this 
approach  was  considered  unrealistic.  The  next  approach 
Involved  reestablishing  the  coordinate  frame  origin  at  each 
control  update  period  to  the  estimated  beam  or  estimated 
target  location.  This  method  was  actually  implemented, 
first  based  on  estimated  target  position  and  then  with 
estimated  beam  position.  The  results  were  virtually 
identical.  Not  only  did  this  approach  maintain  the  beam 
over  the  detector,  it  had  a  marked  impact  on  the  stability 
of  the  closed-loop  system.  Figure  5.12  depicts  the  mean 
error  +/-  1  standard  deviation  envelope  of  the  adaptive 
state  estimate.  However,  it  is  also  representative  of  the 
performance  of  the  first  and  second  Meer  filters  with  beam 
time  constants  of  20  and  30  seconds  respectively.  Figure 
5.13  depicts  the  performance  of  a  Meer  filter  with  a  beam 
time  constant  of  10  seconds.  The  notable  difference  in  this 
filter's  performance  is  the  steady  state  bias  error.  The 
wider  range  of  filter  beam  time  constants  and  a  true  beam 
time  constant  of  25  seconds  was  picked  to  emphasize  the 
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robustness  of  this  approach,  in  Johnson's  research 
[9:119-123]  a  filter-assumed  beam  time  constant  of  20.0 
seconds  and  a  true  beam  time  constant  of  19.2  seconds  or  21 
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seconds  resulted  in  an  unstable  beam  filter.  What  is 
important  in  Figure  5.12  is  that  the  adaptive  mean  is  near 
zero  and  the  filter  computed  error  standard  deviation  (as 
indicated  by  the  two  dashed  lines  that  straddle  zero)  is 
close  to  the  actual  error  standard  deviation.  Note  that  the 
MMAE  filtered  out  the  mean  error  of  the  third  Meer  filter 
(through  its  probability  calculations).  At  the  time  of  the 
actual  runs  for  this  investigation,  the  software  had  not 
been  modified  to  output  the  individual  probabilities  nor  the 
adaptive  beam  estimate.  However,  from  the  results  of  the 
previous  sections,  we  may  surmise  that  the  MMAE  estimate  was 
dominated  by  the  estimate  of  the  Meer  filter  with  the 
largest  beam  time  constant.  Important  in  Figure  5.13  is  the 
slope  of  the  mean  error  in  steady  state  remains  flat  and  the 
error  standard  deviation  is  well  represented  by  the  filter 
computed  value,  despite  the  presence  of  the  steady  state 
bias  error.  For  none  of  these  filters  is  there  an 
Indication  of  instability.  Significant  in  these  figures  is 
the  large  initial  error  and  the  long  transient  before 
reaching  steady  state.  For  this  run,  as  in  all  the  runs, 
the  initial  target  and  beam  positions  were  set  to  zero  in 
both  the  truth  and  filter  models.  The  initial  target 
velocity  state  was  5  cm/sec  and  initial  acceleration  state 
was  also  zero.  The  initial  error  is  due  to  the  filter's 
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lack  o£  information  during  startup.  When  the  initial 
sequence  of  measurements  includes  noise  measurements,  the 
initial  estimate  of  the  beam's  position  is  very  poor.  Since 
all  events  now  occur  on  the  array,  all  measurements  are 
dewelghted  through  Equations  (2-7),  (2-39)  and  (2-46). 

Thus,  it  requires  additional  signal  measurements  for  the 
filter  to  update  toward  the  beam  location.  It  should  be 
noted  that  these  results  were  obtained  by  supplying  the  MMAE 
Meer  estimate  to  the  MMAC  developed  by  Johnson  19)  with  no 
retuning  or  adjustments.  Figure  5.14  indicates  how  well 
this  controller  performed  in  tracking  the  beam  onto  the 
target.  The  upper  line  is  the  true  RMS  error  in  tracking 
error,  whereas  the  lower  line  is  the  filter-computed  RMS 
tracking  error.  Significant  in  this  figure  is  that  the  true 
RMS  tracking  error  is  flat  in  steady  state  and  the  RMS  error 
of  the  filter  is  close  to  this  value.  The  large  initial 
transient  is  due  the  initial  transient  in  the  beam  filter 
estimate . 

Although  this  approach  appeared  to  return  robustness  to 
the  closed-loop  system,  the  author  discontinued  further 
development  along  this  path.  This  was  done  because 
reinitializing  the  coordinate  frame  or  the  use  of  relative 
states  versus  absolute  states  in  the  real  world  represents 
rotating  the  entire  NPB  device  toward  the  target.  In  a 
simulation  this  can  be  done  easily,  instantaneously  and 
without  error.  In  the  real  world  this  may  be  impractical 
for  tracking  a  highly  dynamic  target  due  to  the  mass  of  the 
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NPB  device.  Thus,  this  approach  may  be  more  useful  in 
providing  slow  coarse  pointing  adjustment.  The  third 
approach  investigated  a  more  practical  method  of  tracking  a 
highly  dynamic  target.  Instead  of  moving  the  entire  NPB 
device,  the  detector  is  moved.  This  would  then  provide 
vernier  or  fine  tracking  of  the  target.  For  simplicity  (as 
this  is  a  feasiblity  study),  the  detector  is  repositioned 
instantaneously  at  the  time  of  controller  calculations.  The 
new  position  of  the  detector  is  based  on  the  estimated  beam 
location  halfway  through  the  next  controller  update  cycle. 
This  was  done  in  an  effort  to  ensure  that  the  beam  would  be 
over  the  detector  so  all  the  signal  events  would  originate 
on  the  array.  As  discussed  in  Section  5.6.1  any  signal 
events  "occurring"  off  the  array  result  in  the  filter  being 
artifically  informed  that  the  event  was  a  signal.  In 
essence,  we  are  no  longer  telling  the  filter  which  events 
are  due  to  signal  and  which  are  due  to  noise. 

This  approach  does  not  provide  the  same  stability  found 
with  the  first  approach.  Figures  5.15  and  5.16  indicate  the 
performance  of  the  MMAE  Meer  estimate  and  the  best  Meer 
filter  performance  which  was  provided  by  the  Meer  filter 
with  the  largest  beam  time  constant.  In  this  run,  the  range 
of  beam  time  constants  were  set  at  20,  20.4,  and  19.6 
seconds  with  the  true  beam  time  set  to  20.2  seconds  to 
minimize  the  growth  of  the  actual  error  standard  deviation. 
This  was  necessary  so  the  entire  error  time  history  could  be 
shown  on  one  figure  without  deemphasizing  fine  detail,  as 


would  have  occurred  had  the  error  standard  deviation  been 
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large.  Important  in  these  figures  is  the  poor  filter 
computed  error  standard  deviation,  which  is  much  less  than 
the  actual  error  standard  deviation.  In  addition,  the  mean 
error  does  not  remain  near  zero  for  the  elemental  Meer 


filter  estimate,  nor  the  adaptive  estimate.  Most  important 
is  that  the  actual  error  standard  deviation  is  increasing 
with  time.  Both  the  drift  of  the  mean  error  from  zero  and 
the  growing  error  standard  deviation  indicate  a  divergent 
filter.  The  performance  is  also  depicted  in  Table  5.8. 

Note  that  the  square  root  of  the  filter  predicted  variance 
is  again  near  the  square  root  of  the  steady  state  variance 
of  the  beam  calculated  through  Equation  (5-1).  Thus,  the 
elemental  Meer  filters  are  not  making  effective  use  of  the 
incoming  measurements.  Despite  this,  the  MMAE  Meer  estimate 
was  able  to  combine  the  poor  elemental  estimates  into  a 
marginally  better  estimate.  As  with  the  previous  study. 


Table  5.8  Filter  Performance  with  Detector  Following  Beam 


R  =  . 

True 

2 

5  cm  g 

beam  time 

=  .2  cm/sec 

constant  = 

1/2  Depth  =  1 

20.2  seconds 

beam 

avg  mean 

square  root  of 

t  ime 

constant 

error 

filter  variance 

V^SS 

( sec ) 

(cm) 

(cm) 

Adaptive 

.148876 

.555665 

Meer 11 

20.0 

.277048 

.665026 

.632 

Meer#2 

20.4 

-.170792 

.673304 

.639 

Meer #3 

19.6 

.752915 

.656665 

.626 

20.2 

.636 
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this  run  was  actually  conducted  prior  to  the  modifications 
to  the  code  necessary  to  allow  the  elemental  probability 
histories  to  be  sampled  and  output  so  their  ensemble 
averages  could  be  calculated.  We  can  see  from  Figure  5.15 
that  the  MMAE  estimate  is  diverging,  just  as  the  elemental 
filters  are  diverging.  There  are  two  reasons  for  the 
elemental  filter  divergence.  The  first  reason  is  the 
deweighting  of  the  measurements  by  the  Meer  filter  as  has 
already  be  discussed.  However,  as  we  previously  noted, 
Johnson  [9]  experienced  divergence  even  when  the  Meer  filter 
was  fully  weighting  signal  measurements  and  providing 
neglible  weight  to  noise  measurements.  The  second  reason  is 
due  to  the  physical  nature  of  the  problem  of  tracking  an 
object  that  can  move  to  infinity.  Recall  from  Section  3.1.1 
that  our  target  model  involves  a  double  integration  process. 
To  track  this  target,  the  controller  must  force  the  beam 
state,  as  modelled  in  this  simulation,  to  infinity.  As  the 
beam  state  grows  large,  the  propagation  error  to  the  next 
beam  state  due  to  a  mismodelling  of  the  beam  time  constant 
also  grows  large  as  depicted  in  Figure  5.17.  In  this  figure 
the  upward  ramping  line  is  the  propagation  error  committed 
by  a  filter  with  a  time  constant  of  20.8  seconds.  The 
downward  ramping  line  is  the  propagation  error  committed  by 
a  filter  with  a  time  constant  of  19.2  seconds.  The 
propagation  equation  used  is  Equation  (2-17)  with  a 
propagation  time  interval  of  1  second  and  a  control  input  of 
10.  The  control  input  chosen  represents  a  nominal  value 


B 


necessary  to  track  the  target.  Although  a  control  input  can 
be  applied  to  eliminate  the  error  (a  control  input  equal  to 


the  beam  state)  this  would  not  in  general  track  the  target 
(unless  the  target  position  is  reset  to  zero  as  in  the  first 
approach).  As  seen  in  Figure  5.17,  the  propagation  error 


eventually  exceeds  the  beamwidth  and  so  signal  measurements 
become  deweighted  just  as  noise.  Therefore,  it  would  seem 
that  the  ultimate  solution  to  the  stability  problem  involves 
rotating  the  NPB  device  or  using  relative  states  versus 
absolute  states.  However,  it  may  still  be  possible  to 
extend  operations  beyond  propagation  errors  that  exceed  the 
beamwidth  using  the  MMAE  techniques.  However,  the  number  of 
elemental  filters  would  be  dependent  on  the  beam  position. 
The  concept  is  to  have  sufficient  elemental  filters  such 
that  the  propagation  error  between  adjacent  filters  (in 
parameter  space)  does  not  exceed  the  beamwidth.  The 
elemental  filters  near  the  true  beam  time  constant  would 
make  effective  use  of  the  incoming  measurements  and  remain 
stable.  The  other  filters  would  eventually  become 
divergent.  However,  the  divergent  filters  could  be  reset 
based  on  the  adaptive  estimate  of  the  nondivergent  filters. 
To  do  this,  we  would  need  to  detect  when  an  elemental  is 
divergent.  Figure  5.18  depicts  the  mean  residual 
performance  of  the  Meer  filter  with  a  beam  time  constant  of 
19.6  seconds.  Significant  here  is  positive  growth  in  the 
residuals,  indicating  the  negative  propagation  error  in  this 


filter'3  estimate  of  the  true  beam  position,  as  indicated  by 


Figure  5.17.  Thus,  through  residual  monitoring  we  may  be 
able  to  detect  the  divergence  of  the  filter. 


5.6.2  Divergence  Detection  and  Filter  Reset.  As  the 
elemental  Meer  filters  would  surely  be  diverging,  it  was 
necessary  to  design  a  method  of  determining  divergence  and 
then  to  reset  the  filters.  The  author  chose  to  detect 
divergence  through  residual  monitoring.  As  a  filter 
diverges,  the  mean  of  the  residuals  will  increase,  as  seen 
in  Figure  5.18,  or  their  variances  will  increase,  whereas 
the  filter  expectation  of  the  residual  squared  values  given 
by  (scalar  measurement)  [16:2291, 


E{r2(ti)}  *  H(t1)P(ti)H(t1)  +  (5-2) 

would  remain  at  a  steady  state  value  corresponding  to  the 
steady  state  value  of  the  variance  P(tj).  What  is  needed  is 
a  method  to  determine  when  a  filter  has  diverged,  so  that  it 
can  be  reset.  Since  we  have  assumed  the  residual  sequence 
is  a  zero-mean  Gaussian  stochastic  process,  we  will  detect 
divergence  with  the  likelihood  function  [16:230], 

LM(t.)  =  c  ( t ,  )  -  1/2  E  r  2  ( t  . )  /  cr2(t.)  (5-3) 

N  1  1  j=i -N+l  j  j 

where  N  is  the  number  of  residuals  being  processed,  c  is  an 
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arbitrary,  slowly  varying  value  and 

l/or2(tj)  =  IH(t  +  RUj)]"1  (5-4) 

For  this  design  c  Is  set  to  zero.  The  output  of  the 
likelihood  function  Is  checked  against  a  user  selected 
threshold  that  Indicates  divergence.  For  this  study  the 
design  consisted  of  individual  divergence  detectors 
(likelihood  functions)  for  each  filter  and  reset  a  divergent 
filter  individually  when  a  threshold  of  9  (three  standard 
deviations  squared)  was  exceeded.  As  a  filter  within  the 
bank  diverges,  its  probability  should  drop  relative  to  the 
nondivergent  filters.  The  impact  of  this  filter  on  the 
adaptive  state  estimate  then  becomes  small.  Thus,  for  this 
study,  the  adaptive  state  estimate  is  used  to  reset  a 
divergent  filter.  Properly,  the  divergence  check  should 
occur  prior  to  the  calculation  of  the  adaptive  state 
estimate  so  the  divergent  filters  can  be  removed  from  the 
estimation  process.  However,  this  correction  was  not 
incorporated  into  this  study.  This  design  used  a  value  of 
N  =  5  in  Equation  (5-3).  This  value  was  chosen  by 
determining  the  approximate  time  that  the  state  estimate  of 
an  elemental  filter  with  beam  time  constant  of  19.2  seconds 
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differed  from  the  true  state  by  three  beam  halfwidths  when 
the  true  beam  time  constant  is  20.4  seconds.  The  value  of  N 
was  then  set  so  resetting  of  the  filter  occurred  shortly 
after  that  time.  Although  crude,  the  author  deemed  this 
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adequate  £or  a  signal  to  noise  ratio  of  ten.  For  lower 
signal  to  noise  ratios,  a  larger  value  of  N  would  be  more 
appropriate.  A  larger  value  of  N  would  reduce  the  chance  of 
improperly  resetting  the  filters  based  on  consecutive  noise 
events.  Recall  that  noise  events  are  equally  likely  over 
the  entire  array  rather  than  clustered  near  a  mean 
regardless  of  the  performance  of  the  filter.  Recall  also 
that  signal  events  can  also  have  large  residuals  for  a 
nondivergent  filter,  although  the  arrival  rate  is  low  for 
such  events.  Unfortunately,  a  larger  value  for  N  would 
delay  the  detection  and  reset  of  a  divergent  filter.  Thus, 
the  adaptive  estimate  may  suffer.  Care  must  also  be  taken 
in  selecting  a  threshold  for  declaring  divergence  in  a  noisy 
environment.  Since  the  noise  rate  is  constant  over  the 
array,  noise  events  at  the  array  edges  have  the  same 
probability  of  occurring  as  noise  events  near  the  beam. 

Also,  we  have  defined  the  minimum  appropriate  detector 
length  in  Section  2.2.1  as  12  beamwidths.  Thus,  most  of  the 
noise  events  are  considered  to  have  large  residuals  even  for 
a  nondivergent  filter.  Therefore,  a  low  threshold  can  be 
exceed  with  a  few  noise  events  despite  a  large  value  of  N. 
This  drives  a  desire  to  minimize  the  array  length  to  limit 
the  size  of  residuals  from  noise  events,  as  well  as,  the 
rate  of  noise  events. 

5.6.3  MMAE  Meer  Baseline.  With  the  Meer  probability 
bound  set  to  .0005  as  mentioned  in  Section  5.2  and  the 


filter  divergence  detection  and  reset  processes  in  place, 
three  runs  were  made  to  establish  a  baseline  of  Meer  filter 
state  estimation  and  parameter  estimation  performance  in  a 
closed  loop.  The  controller  for  this  closed  loop  baseline 
was  designed  with  the  weighting  matrices  Xu  =  100  (target 
tracking  error  weight),  X  =  10  (pseudointegral  state 
weighting),  and  U  =  1  (control  weighting)  used  by  Johnson 
19].  Since  the  desire  of  this  section  is  to  investigate  the 
MMAE  Meer  operation  and  not  the  controller  operation,  the 
controller  was  provided  the  true  beam  and  target  states  for 
these  runs.  The  Meer  filters  1,  2,  and  3  beam  time 
constants  were  20,  20.8,  and  19.2.  The  first  run  had  a 
sinusoidal  beam  time  constant  with  an  initial  value  of  20 
seconds  and  completed  one  cycle  in  100  seconds  as  in  Section 
5.4.  Since  preliminary  runs  indicated  that  the  adaptive 
parameter  estimation  process  was  functioning  properly,  the 
propagation  of  the  adaptive  state  estimate  to  measurement 
times  used  the  adaptive  parameter  estimate  provided  by  the 
filter.  Figure  5.19  depicts  the  mean  error  of  the  adaptive 
estimate.  Figure  5.20  depicts  the  response  of  the  Meer 
filter  based  on  a  beam  time  constant  of  20  seconds.  It  is 
representive  of  the  performance  of  the  elemental  filters. 
Both  the  MMAE  Meer  filter  and  the  elemental  Meer  filters 
appear  to  be  underestimating  their  own  errors.  Also  the 
actual  error  standard  deviations  are  increasing.  This 
Indicates  the  filters  are  continuing  to  diverge.  As  we  have 
provided  the  true  states  to  the  controller,  we  have  isolated 


closed  loop  beam  errors  to  the  beam  filter.  Note  that  the 
parameter  estimation  performance  seen  in  Figure  5.21  has 
improved  from  the  open  loop  case  as  indicated  by  the  general 
result  of  Figure  5.3.  The  dots  on  Figure  5.21  serve  as 
reference  points  of  the  correct  beam  time  constant.  As  can 
be  seen,  the  adaptive  beam  time  constant  is  a  reasonably 
good  estimate  of  the  true  beam  time  constant,  with  some  lag. 
The  adaptive  state  estimate,  in  this  case,  is  far  better 
than  any  of  the  elementals  as  one  might  expect. 

The  second  baseline  run  used  a  steady  true  beam  time 
constant  of  20.4  seconds.  Figures  5.22  and  5.23  depict  the 
MMAE  performance  and  a  representative  performance  of  the 
elemental  estimators.  Again  divergence  in  both  figures  is 
noted.  Figure  5.24  shows  the  adaptive  parameter  estimate. 
Although  the  mean  adaptive  estimate  is  close  to  the  correct 
value  of  20.4  seconds,  it  is  subject  to  variation.  Again, 
this  is  reasonable  parameter  estimation  performance.  The 
mean  probability  time  histories  are  shown  in  Figures  5.25, 
5.26,  and  5.27.  Significant  in  these  figures  is  that  the 
mean  probability  is  lowest  for  the  Meer  filter  with  beam 
time  constant  that  is  the  farthest  from  the  true  beam  time 
constant.  Overall  the  parameter  estimation  performance  of 
the  MMAE  Meer  filter  for  a  controlled  beam  is  good. 

The  third  run  used  the  nominal  beam  time  constant  of  20 
seconds.  This  value  was  chosen  to  examine  the  effect  of 
incorporating  the  MMAE  Meer  filter  structure  with  divergence 
detection  and  filter  reset  processes  when  the  true  beam  time 
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constant  matches  one  of  the  filter -assumed  beam  time 
constants.  Figure  5.28  shows  the  performance  of  the  MMAE. 

As  before,  the  actual  standard  deviation  exceeds  the  filter 
computed  standard  deviation.  In  addition,  the  error 
standard  deviation  of  the  adaptive  estimate  is  growing. 

This  is  rather  surprising  given  that  one  of  the  Meer  filter 
assumed  beam  time  constant  matches  the  true  value.  However, 
Figure  5.29  indicates  that  even  the  filter  that  matches  the 
beam  time  constant  is  underestimating  its  error.  It  might 
appear  that  the  filter  requires  some  tuning.  However,  it 
should  be  noted  that  the  filter  dynamics  model  matches  the 
true  dynamics  model  and  the  filter  parameter  values  R  and  g 
match  those  of  the  truth  model.  Figure  5.30  is  the  RMS  beam 
error  time  histories  (true  and  filter-computed)  of  the  Meer 
filter  with  the  correct  beam  time  constant.  The  upper  solid 
line  is  the  true  beam  error  committed  by  this  filter.  The 
lower  broken  line  is  filter-computed  one  standard  deviation. 
This  figure  clearly  indicates  that  the  filter  is  under¬ 
estimating  its  error.  Since  the  actual  error  of  this  filter 
in  growing  with  time  the  filter  is  divergent.  As  we  are 
supplying  the  controller  with  the  true  states,  we  have  again 
isolated  the  error  to  the  beam  filter.  As  discussed 
previously,  the  filter  once  again  appears  to  be  under¬ 
weighting  its  measurements,  resulting  in  divergent  behavior 
in  state  estimation.  Examining  the  mean  adaptive  beam  time 
constant  in  Figure  5.31,  we  note  that  the  adaptive  estimate 
does  remain  near  the  true  value  of  20.  Also,  in  Figure  5.32 
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we  see  the  highest  probability  is  correctly  associated  with 
the  filter  with  beam  time  constant  of  20.  Figures  5.33  and 
5.34  indicate  that  both  of  the  other  filters  have  equally 
low  probability  weights.  Again,  reasonably  good  per¬ 
formance.  However,  we  might  have  expected  the  probability 
weights  of  the  filters  with  the  incorrect  beam  time  constant 
to  be  driven  to  the  lower  bound  of  0.0005.  However,  the 
occurrence  of  a  noise  event  can  result  in  shifts  of  the 
probability  weights  for  the  individual  Meer  filters. 

To  investigate  the  degree  that  the  Meer  filter  is 
deweighting  the  measurements,  a  Meer  filter  was  operated 
with  the  correct  beam  time  constant,  initial  conditions 
matching  the  truth  model,  and  no  measurements.  As  would  be 
expected,  the  actual  error  standard  deviation  envelope  in 
Figure  5.35  closely  matches  the  filter  predicted  standard 
deviation,  which  is  the  square  root  of  steady  state  error 
variance  given  by  Equation  (5-1).  Figure  5.36  is  the 
performance  of  a  single  Meer  filter  with  measurements 
provided.  What  is  significant  between  these  two  figures  is 
the  similarity  between  actual  error  standard  deviations  and 
the  steady  state  variance  of  the  beam  estimate.  They  are 
nearly  the  same  in  steady  state.  This,  once  again, 
indicates  that  the  Meer  filter  is  improperly  deweighting 
measurements  and  is  relying  on  its  internal  dynamics  model. 
The  reason  for  this  is  seen  in  Fiqure  2.4,  the  slopes  of  the 
signal  rate  function  are  rather  steep.  Thus,  a  small 
variation  of  the  filter-estimated  beam  position  from  the 


actual  beam  position  can  result  in  a  low  estimate  of  the 
signal  rate  and  the  deweighting  of  elemental  Snyder-Fishman 
filters  that  assumed  the  event  to  be  a  signal.  The  reason 
we  see  this  in  the  closed  loop  case  and  not  the  open  loop 
case  is  the  fact  that  in  open  loop  the  truth  model  and  the 
filter  models  will  propagate  toward  each  other  if  the  sign 
of  the  true  beam  position  is  opposite  that  of  the  filter 
assumed  position  (all  propagate  toward  zero),  whereas  in 
closed  loop,  the  sign  of  the  true  beam  position  and  the 
f ilter-assummed  beam  position  will  generally  be  the  same  as 
the  beam  is  driven  toward  the  target  and  generally  away  from 
zero.  Again  this  is  a  result  of  operating  with  absolute 
states  rather  than  relative  states.  Note  In  Figure  5.36, 
that  the  error  standard  deviation  envelope  is  not  growing. 
This  demonstrates  that  a  Meer  Filter  operating  outside  the 
MMAE  Meer  structure,  with  the  correct  beam  time  constant, 
and  starting  at  the  same  initial  conditions  as  the  truth 
model.  Is  at  least  stable.  This  is  what  we  would  expect. 
When  compared  to  Figure  5.29  we  see  a  possible  indication 
that  noise  events  are  triggering  the  unnecessary  reset  of 
this  filter  while  In  the  MMAE  Meer  structure. 

5.6.4  Controller  Baseline.  To  establish  the  controller 
baselines,  the  performance  of  the  controller  with  full  state 
feedback  that  was  previously  discussed  and  the  performance 
of  the  same  controller  with  estimates  provided  by  a  single 
Meer  and  a  single  Kalman  filter  are  depicted  in  Figures  5.37 


and  5.38.  The  Meer  filter  and  Kalman  filter  both  use  the 
true  parameter  values.  As  expected,  the  true  RMS  error 
(Figure  5.37  lower  line)  for  full-state  feedback  is  very 
small.  The  upper  line  is  the  RMS  error  as  calculated  from 
the  error  between  the  MMAE  Meer  estimate  and  the  target 
Kalman  filter  based  on  the  true  dynamics  driving  noise 
strength.  In  Figure  5.38  we  see  the  effects  of  using  filter 
estimates  on  tracking  performance.  Now  the  upper  line  is 
the  true  RMS  error  and  the  lower  line  is  the  RMS  error  based 
on  filter  estimates.  As  one  would  expect,  the  RMS  error 
increases  when  full-state  estimation  is  replaced  with  filter 
state  estimates.  Table  5.9  contains  the  RMS  tracking  errors 
for  these  two  baselines.  The  RMS  tracking  error  of  the  nine 
element  MMAC  is  depicted  in  Figure  5.39.  Here,  the  upper 
line  is  the  actual  tracking  error  whereas  the  lower  line  is 
the  RMS  tracking  error  as  calculated  from  the  MMAE  Meer 
estimate  and  the  Kalman  filter  estimate  based  on  the  true 
dynamics  driving  noise  strength  of  0.01.  Note  that  the  RMS 
error  is  increasing  with  time,  indicating  divergence.  Such 


Table  5.9  Controller  Baselines 


2  1/2 

R  =  . 5  cm  g  =  .2  cm/sec  Depth  =  1 

True  beam  time  constant  =  20  seconds 


controller 

type 

Full-state  Feedback 
Filter  Estimates 


RMS  tracking  error 
(cm) 

.212571 

1.60719 


poor  performance  may  be  due  to  the  Kalman  filter  lower  bound 
of  .1.  This  bound  results  In  poor  elemental  controllers 
receiving  a  greater  probability  weight  than  would  normally 


occur . 

5.7  Implicit  Model  Following  Controller  Performance 

Programming  errors  coupled  with  time  constraints 
prevented  analysis  of  the  Implicit  Model  Following 
controller  design. 

5.8  Summary 

This  chapter  has  presented  the  results  of  the  analysis 
described  in  Chapter  4.  Performance  of  the  MMAE  Meer  filter 
with  a  wide  variety  of  operating  conditions  was  presented. 
While  in  open  loop,  the  MMAE  Meer  filter  is  able  to  estimate 
the  beam  state  with  little  difficulty.  However,  it  is 
unable  to  estimate  the  beam  time  constant  in  open  loop 
operation.  The  adaptive  beam  time  contant  estimate  tends 
toward  the  largest  beam  time  constant.  This  is  attributed 
to  the  superior  residual  performance  of  a  Meer  filter  with  a 
large  beam  time  constant.  Open-loop  performance  of  the  MMAE 
Meer  filter  prompted  an  investigation  of  the  closed-loop 
performance.  This  resulted  in  the  examination  of  techniques 
to  keep  the  beam  over  the  array.  One  approach  was  to  rotate 
the  NPB  toward  the  target.  The  results  of  this  approach 
indicates  that  the  measurement  are  being  deweighted  by  the 


Meer  filter.  The  effect  of  the  deweighting  is  an  extended 
initial  transient.  However,  steady  state  performance  was 
good.  In  another  approach,  the  detector  is  periodically 
repositioned  under  the  beam.  Although  no  long  transient 
response  is  noted,  there  is  filter  divergence.  To  counter 
this  characteristic,  a  likelihood  function  is  used  to  detect 
divergence.  When  it  diverges  a  Meer  filter  is  reset  to  the 
MMAE  Meer  filter  state  estimate.  Under  these  condition  the 
state  estimation  performance  is  still  divergent.  However, 
there  is  improvement  in  the  parameter  estimation  performance 
of  the  MMAE  Meer  filter.  When  this  parameter  is  used  to 
propagate  the  adaptive  state  estimate,  the  performance  of 
the  MMAE  state  estimation  is  superior  to  the  elemental 
estimates.  However,  analysis  indicates  that  again  the 
measurements  are  being  dewelghted  in  the  elemental  Meer 
filters.  Throughout  the  analyses,  there  is  seen  a  dramatic 
difference  in  state  estimation  performance  and  parameter 
estimation  performance  between  the  open  loop  case  and  the 
closed  loop  case.  This  difference  is  also  seen  in  the 
different  methods  used  to  track  the  target.  This  is 
attributed  to  the  difference  between  using  relative  states 
and  absolute  states. 

Two  performance  baselines  are  established  for  a 
controller  designed  around  the  the  weighting  matrices 
Xu  *  100,  X^  =  10,  and  U  -  1  as  used  by  Johnson.  The 
first  baseline  is  based  on  this  controller  receiving 
full-state  feedback.  The  second  baseline  is  based  on  this 


same  controller  receiving  filter  estimates  from  a  single 
Meer  filter  and  a  single  Kalman  filter  based  on  the  true 
parameter  values.  The  RMS  tracking  error  for  a  nine-element 
MMAC  is  presented.  This  controller  appeared  to  be  more 
unstable  than  a  three-element  MMAC  receiving  adaptive  beam 
state  estimates  from  the  MMAE  Meer. 


VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


6 . 1  Conclusions 

Eight  analyses  of  MMAE  Meer  filter  and  a  nine-element 
MMAC  controller  were  completed.  The  MMAE  Meer  filter  for 
this  study  was  composed  of  3  elemental  Meer  filters  based  on 
separate  assumed  beam  time  constants.  The  first  five 
analyses  studied  MMAE  Meer  filter  performance  with  an 
uncontrolled  beam.  The  last  three  analyses  examined  the 
MMAE  Meer  filter  performance  with  a  controlled  beam.  As 
part  of  the  closed  loop  analyses,  baselines  with  the 
three-element  MMAC  used  by  Johnson  [9]  were  established. 
Finally,  the  RMS  tracking  performance  of  the  nine-element 
MMAC  was  presented.  The  first  analysis  indicates  that 
there  was  little  change  in  state  estimation  performance  for 
a  small  mismodelllng  of  the  beam  time  constant  for  an 
uncontrolled  beam.  Results  of  the  second  analysis  show  that 
the  MMAE  Meer  filter  was  unable  to  provide  an  accurate 
parameter  estimate  for  an  uncontrolled  beam.  Although  there 
was  no  major  computational  savings  realized  in  this 
simulation  for  operating  at  a  Meer  filter  depth  of  one  due 
to  the  computational  load  of  other  portions  of  the  software, 
the  third  analysis  demonstrated  that  a  lower  filter  depth 
was  desirable  based  on  the  RMS  error  performance  for  beam 
state  estimation.  As  in  the  previous  analyses,  the  fourth 


analysis  indicates  that  parameter  estimation  performance  for 
an  uncontrolled  beam  was  poor  when  the  beam  time  constant 
was  linearly  or  sinusoidually  varying  in  time.  The  fifth 
analysis  showed  that  state  estimation  performance  was 
sensitive  to  very  large  mismodelling  of  the  beam  time 
constant  for  an  uncontrolled  beam.  The  results  of  the  first 
five  analyses  indicate  that  an  MMAE  Meer  filter  was 
unnecessary  for  tracking  an  uncontrolled  NPB,  i.e.  that  a 
conventional  Meer  filter  is  sufficient. 

However,  the  results  of  the  last  three  analyses 
demonstrate  that  state  estimation  of  the  controlled  beam  is 
very  sensitive  to  variations  in  the  filter  assumed  beam  time 
constant.  In  this  case,  the  MMAE  Meer  filter  is  effective 
In  providing  a  coarse  estimate  of  the  beam  time  constant, 
but  there  appears  to  be  divergence  in  the  beam  state 
estimates.  The  analysis  of  this  divergence  indicates  that 
the  Meer  filter  is  underweighting  the  measurement  data. 

There  was  no  evidence  of  a  drastic  loss  of  control;  however, 
growth  in  RMS  tracking  error  is  evident.  The  final  analysis 
demonstrated  that  the  controller  is  sensitive  to  poor  state 
estimates  provided  by  the  filters.  The  first  baseline 
controller  performance  indicated  that  the  controller  design 
performed  well  with  full-state  feedback.  Incorporation  of 
filter  estimates  from  filters  based  on  the  true  parameters 
in  the  second  baseline  resulted  in  a  substantial  Increase  in 
rms  tracking  errors.  However,  the  RMS  tracking  errors  did 
not  Increase  with  time,  indicating  at  least  stability  at 


nominal  conditions. 

Analysis  of  the  closed  and  open  loop  performances 
indicate  there  is  a  substantial  improvement  in  state 
estimation  performance  when  state  values  are  maintained  near 
zero  (using  relative  states  versus  absolute  states).  The 
parameter  estimation  performance  improved  in  the  closed  loop 
case  as  the  propagation  differences  between  the  elemental 
Meer  filter  led  to  significant  variations  in  residual 
performance.  There  is  also  an  indication  that  divergence 
detection  based  on  a  likelihood  function  with  a  threshold  of 
nine  and  N  =  5  results  in  Improper  resetting  of  nondivergent 
filters . 

No  performance  data  was  collected  on  a  controller  based 
on  implicit  model  following  techniques.  However*  there  are 
conclusions  that  can  be  drawn  from  the  development  process. 
Implicit  model  following  is  not  well-suited  for  problems 
with  a  first  order  Gauss-Markov  plant,  with  one  control 
input,  and  a  single  output.  This  reduces  the  model  that  can 
be  followed  to  a  single  dimension  (see  Section  3. 2. 3. 2). 

This  does  not  offer  a  great  deal  of  flexibility  in  closed 
loop  pole  selection.  It  must  be  noted  that  other  techniques 
are  also  limited  by  the  same  constraint.  Given  a  desirable 
closed  loop  pole  based  on  design  specifications  for 
performance  and  robustness,  this  method  does  provide  a 
straightforward  approach  to  matching  that  pole. 


6 . 2  Recommendations 

These  recommendations  are  divided  into  three  parts.  The 
first  part  concerns  continued  research  with  the  beam 
estimator.  The  second  concerns  controller  development.  The 
last  involves  other  avenues  of  research  and  improvement  of 
the  research  tools. 

Of  foremost  importance  to  the  continued  use  of  the  Meer 
filter  in  beam  control  is  the  demonstration  that  the  Meer 
filter  can  be  tuned  for  closed  loop  operation.  As  indicated 
in  Chapter  V,  the  Meer  filter  apparently  well  tuned  for  open 
loop  operations,  appears  to  be  poorly  tuned  for  closed  loop 
operation.  Once  tuning  capability  is  demonstrated,  it  will 
be  necessary  to  repeat  sensitivity  and  robustness  studies 
for  variations  in  the  beam  parameters.  The  first  step  in 
attacking  this  tuning  problem  should  be  a  sensitivity 
analysis  of  the  Meer  filter  with  respect  to  the  beam 
dispersion  R,  with  the  detector  properly  located  underneath 
the  beam,  during  closed  loop  operations.  Either  of  two 
methods  can  be  considered  for  maintaining  the  beam  over  the 
detector.  If  rotation  of  the  coordinate  frame  (using 
relative  variables)  is  used,  then  the  emphasis  will  be  on 
tuning  for  transient  response.  If  the  detector  is  swept 
under  the  array  (using  absolute  variables),  then  the 
emphasis  will  be  on  the  steady  state  behavior.  The  goal  of 
either  approach  will  be  to  allow  the  Meer  filter  to  make 
more  effective  use  of  measurement.  However,  if  the  second 


approach  is  selected.  It  is  recommended  that  the  detector  be 
continuously  repositioned  under  the  array  instead  of 
periodically  as  done  in  this  study.  Properly  locating  the 
detector  may  necessitate  modelling  the  detector  motion  and 
providing  appropriate  control  to  ensure  all  the  signal 
events  occur  on  the  detector  array.  The  particular  values 
of  R  and  square  root  of  beam  dynamics  driving  noise  g  used 
in  this  study  may  have  contributed  greatly  to  the  apparent 
deweighting  of  measurements.  The  results  of  the  sensitivity 
study  may  indicate  a  more  appropriate  value  or  range  of 
values  of  R  for  subsequent  study.  The  second  step  may 
involve  the  purposeful  detuning  of  the  Meer  filter  at  open 
loop  conditions  to  provide  better  performance  for  closed 
loop  conditions.  Increasing  the  square  root  of  the  dynamics 
driving  noise  g  in  the  filter  is  the  simplest  method.  This 
would  boost  the  covariance  and  provide  additional  gain  to 
the  measurement.  However,  such  a  gain  would  be  applied 
indiscriminately  and  may  result  in  undesirable  performance 
in  noisy  environments. 

Although  the  MMAE  Meer  filter  was  effective  in  providing 
a  coarse  beam  time  constant  estimate  in  the  feedback  mode,  a 
much  finer  estimate  may  be  possible  with  a  finer 
discretization  of  the  beam  time  constant.  This  would 
involve  placing  additional  Meer  filters  within  the  MMAE  Meer 
structure.  This  may  also  reduce  the  phase  lag  in  estimating 
a  sinusoidually  varying  beam  time  constant.  Therefore  it  is 
suggested  that  performance  of  the  MMAE  Meer  with  additional 


elemental  Meer  filters  be  evaluated.  If  absolute  variable 
are  used  then  it  also  recommended  that  the  number  of 
elemental  filters  used  in  the  MMAE  Meer  become  a  function  of 
the  beam  state  estimate.  In  addition,  the  state  estimate 
used  to  reset  divergent  filters  should  be  based  on  the 
nondivergent  filters. 

One  aspect  of  the  MMAE  Meer  filter  that  needs 
Investigation  is  Its  performance  under  very  noisy 
conditions.  The  probability  weights  used  to  determine  the 
adaptive  state  estimate  from  the  individual  Meer  filter 
state  estimates  are  based  on  the  measurement  history  or 
residual  performance.  Since  noise  events  are  a  part  of  that 
residual  performance,  they  may  have  a  substantial  Impact  on 
performance.  When  a  noise  event  occurs  It  Increase  the 
probability  of  poorly  performing  filters.  For  low  noise 
conditions,  the  subsequent  signal  events  would  cause  the 
probabilities  to  fall  for  filters  that  are  not  accurately 
predicting  the  beam  location.  In  a  very  noisy  environment. 
It  is  conceivable  that  a  sequence  of  noise  events  could 
occur  that  would  result  in  a  high  probability  weight  being 
assigned  to  a  filter  that  is  actually  providing  a  poor 
estimate  of  the  beam  location.  Such  a  sequence  of  noise 
events  could  also  result  In  the  resetting  of  filters  that 
are  providing  good  estimates  of  the  beam  position,  but  have 
very  high  residuals  with  respect  to  the  noise  events.  This 
could  result  In  loss  of  beam  tracking. 

Another  aspect  of  the  MMAE  Meer  filter  that  needs 
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investigation  and  perhaps  some  development,  is  the  ability 
of  the  MMAE  Meer  filter  to  provide  its  own  estimates  of  the 
maximum  amplitude  of  the  signal  rate  function,  r,  and  the 
noise  rate.  Realistically,  we  will  only  be  able  to  estimate 
the  actual  r  from  our  assumed  beam  dispersion  and  a  guess 
about  the  number  of  photons  that  would  arrive  over  a  given 
time  interval.  The  Meer  filter  decision  process  may  then  be 
called  upon  to  update  this  estimate  based  on  its 
determination  of  the  number  of  signal  events  that  have 
occurred  over  a  given  period  of  time.  Similarly,  the  Meer 
filter's  determination  of  which  events  were  noise  events 
could  be  used  to  provide  an  estimate  of  the  noise  rate. 
Changes  in  laser  and  detector  efficiencies  provide  the 
physical  motivation  to  seek  a  real-time  estimate  of  r  and 
the  noise  rate  for  use  by  the  Meer  filter.  Therefore,  this 
estimation  capability  needs  to  be  tested  and  perhaps 
developed . 

Another  area  of  estimator  development  that  needs 
consideration  is  the  likelihood  function  for  divergence 
detection.  As  previously  mentioned,  the  occurrence  of 
successive  noise  events  may  incorrectly  indicate  divergence 
of  a  filter  that  is  actual  tracking  the  beam  quite 
precisely.  This  filter  would  then  be  unnecessarily  reset. 
However,  it  rc  .y  be  possible  avoid  the  false  divergence 
indication  by  varying  the  setting  of  the  threshold  and  depth 
of  the  likelihood  function  based  on  the  current  estimate  of 
r  and  the  noise  rate.  The  objective  would  be  to  provide 


adequate  divergence  checking  without  becoming  susceptible  to 
the  noise.  However,  this  must  be  weighed  against  the  impact 
of  adaptively  estimating  to  many  parameter.  This  could 
result  in  identif lability  problems  and  possible 
instabilities . 

Although  little  was  accomplished  in  analysis  of  an 
Implicit  Model  Following  controller  design  in  this  study, 
some  recommendations  are  nevertheless  made  for  further 
study.  To  satisfy  sufficient  conditions  for  a  solution  to 
the  Riccati  equation  for  controller  gain  synthesis,  it  was 
nescessary  to  shift  two  poles  of  the  target  model  used  for 
this  synthesis.  It  is  possible  to  avoid  this  modification 
to  the  target  model  used  for  gain  synthesis  by  incorporating 
the  target  model  into  a  Command  Generator  Tracker  [18] 
developed  as  a  precompensator  feeding  into  a  Proportional 
plus  Integral  controller  for  the  beam. 

In  addition,  no  provision  currently  exists  in  the 
simulation  code  to  isolate  the  impact  of  providing  MMAE  Meer 
estimates  and  Kalman  filter  estimates  to  the  PI  controller. 
Apparently,  this  code  was  removed  after  Zicker  completed  his 
study  of  the  PG  controller  with  the  first-order  target  [32]. 
As  targets,  filters,  and  controllers  become  more 
sophisticated,  the  ability  to  isolate  performance  problems 
to  a  specific  part  of  the  system  will  become  increasingly 
important.  It  is  recommended  that  this  capability  be 
returned  to  the  simulation.  This  isolation  process  can  be 
divided  into  four  phases.  The  first  phase  is  full-state 


feedback.  The  second  phase  uses  MMAE  Meer  filter  feedback 
and  target  full-state  feedback.  This  would  allow  isolation 
of  problems  strictly  due  to  the  use  of  MMAE  Meer  filter 
estimates.  It  may  also  provide  insight  into  appropriate 
methods  to  recover  robustness.  Similarily,  phase  three 
would  use  full-state  beam  feedback  and  Kalman  filter 
estimates  to  isolate  and  correct  problems  in  this  area.  The 
final,  fourth  phase  would  use  both  MMAE  Meer  filter  and 
Kalman  filter  estimates  to  determine  what  problems  arise 
from  the  interaction  of  the  estimates  within  the  controller. 
This  can  also  provide  various  baselines  for  controller 
performance . 

As  Johnson  recommended  [9:167],  it  is  also  recommended 
here  that  the  number  of  runs  necessary  to  provide  sufficient 
convergence  of  sample  statistics  to  true  underlying 
statistics  for  controller  evaluation  should  be  rechecked. 

It  is  recommended  that  this  be  accomplished  by  embedding  the 
equations  used  in  SOFEPL  [6:Sec  4]  to  provide  the  Monte 
Carlo  statistics  in  a  user  defined  subroutine  in  SOFE. 
Appropriate  statistics  to  determine  convergence  of  the 
sample  statistics  to  underlying  statistics  could  be 
calculated  at  the  end  of  each  run  using  the  data  from  all 
the  runs  to  that  time  in  the  Monte  Carlo  simulation.  The 
number  of  runs  necessary  for  a  desired  level  of  convergence 
could  then  be  determined. 

The  last  two  recommendations  concern  the  reduction  of 
the  SOFE  computational  load.  With  the  recommendation  to 


incorporate  addition  Meer  filters  into  the  MMAE  Meer 
structure  comes  the  problem  of  increasing  the  computer  time 
required  to  perform  each  analysis.  The  first  recommendation 
is  to  replace  the  differential  equations  used  to  propagate 
the  truth  and  filter  models  with  difference  equations  based 
state  transition  matrices.  The  linear  models  currently  used 
do  not  require  differential  equations  to  provide  precise 
propagation.  If  the  next  study  desires  or  requires  the  use 
of  the  differential  equations,  then  it  is  recommended  that 
the  differential  equation  of  the  dynamics  driving  noise 
strength  be  incorporated  into  the  truth  model  as  suggested 
in  the  SOFE  manual  [24:161.  This  would  allow  the  separate 
noise  injection  period  to  be  eliminated.  Instead,  noise 
would  be  injected  into  the  truth  model  at  the  time  of  Meer 
filter  update  and  Kalman  filter  update  only.  Use  of  the 
differential  equation  would  allow  the  appropriate  discrete 
time  noise  to  be  determined  despite  the  variable  time 
interval.  This  method  is  not  expected  to  reduce  the 
required  computer  time  as  much  as  the  first  method,  but  as 
the  current  noise  Injection  period  is  0.2  seconds,  it  may 
reduce  the  number  of  calls  to  the  integration  routine  by  a 
factor  of  3  by  only  integrating  to  each  Poisson  event  time 
(average  of  one  per  second  in  this  study)  and  to  the  Kalman 
filter  update  time  (once  per  second). 
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